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FOREWORD 

This  report  is  a  review  of  the  fundamental  procedures  developed  for 
solving  practical  USAF  problems  by  means  of  Computational  Fluid  Dynamic 
techniques.  It  represents  the  efforts  of  the  Computational  Aerodynamics 
Group  over  a  several  year  period,  directed  by  Dr.  Wilbur  L.  Hankey  under 
project  2307N603. 
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SECTION  I 

HISTORICAL  BACKGROUND 

First,  a  historical  review  (Reference  1)  of  aerodynamic  theory  shall  be 
accomplished  in  order  to  obtain  an  overall  perspective.  In  1768 
D'Alembert,  an  experimentalist  for  the  King's  Navy,  used  the  potential 
equation  of  fluid  mechanics  and  showed  that,  contrary  to  experimental 
evidence,  zero  force  was  predicted  on  an  arbitrary  body  immersed  in  a 
moving  fluid.  He  stated,  "this  theory  gives  absolutely  zero  resistance: 
a  singular  paradox  I  leave  mathematicians  to  explain."  The  "D'Alembert 
paradox"  lasted  over  a  century  until  1906  when  Joukowski  introduced  the 
concept  of  circulation  (artifically  representing  viscous  effects)  to 
produce  lift.  The  Wright  Brothers  flew  three  years  before  the  Joukowski 
discovery  and  obviously  used  lift.  Prandtl's  boundary  layer  concept  was 
employed  by  Blasius  in  1908  to  analytically  predict  friction  drag. 
Charles  Lindbergh  first  crossed  the  Atlantic  in  1927  and  Chuck  Yeager 
flew  the  Bell  X-l  supersonical ly  in  1947  which  clearly  demonstrates  the 
rapid  development  of  aviation  during  that  period. 

However,  the  first  large-scale  aerodynamic  calculation  was  Kopal's 
(Reference  2)  solution  for  supersonic  flow  over  cones  in  1947.  The  308 
computed  cases  required  three  years  for  a  five-girl  team  to  accomplish 
using  mechanical  adding  machines.  In  1964  Blottner  (Reference  3) 
developed  a  stable  algorithm  to  solve  generalized  boundary  layers  on  an 
electronic  computer.  Two  years  prior  to  this,  John  Glenn  orbited  the 
earth  achieving  a  major  breakthrough  in  our  aviation  history.  In  1971, 
MacCormack  (Reference  4)  published  the  first  significant  numerical 
solution  of  the  compressible  Navier-Stokes  equations. 

By  contrasting  the  flight  achievements  with  theoretical  aerodynamic 
development  two  major  points  can  be  made. 

1)  With  aerodynamic  theory  dramatically  lagging,  how  were  the 
flight  achievements  possible?  The  answer  is  by  means  of  the  wind  tunnel. 
Excellent  wind  tunnel  facilities  provided  the  necessary  data  for  the 
past  advancements  in  aviation. 
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2)  Since  the  Navier-Stokes  equations  have  been  around  since  1827 
why  are  we  only  now  able  to  solve  them?  The  answer  is  that  the  solution 
of  the  Navier-Stokes  equations  required  the  invention  of  a  large-scale 
computer.  The  computer  only  became  "of  age"  during  the  past  decade. 
Prior  to  this  time  only  limited  analytic  calculus  solutions  of 
approximate  equations  could  be  obtained.  Today,  a  "brute  force"  method 
is  used  to  numerically  solve  the  Navier-Stokes  equations  (Reference  5). 
The  method  of  Computational  Fluid  Dynamics  (CFD)  is  based  on  the  fact 
that  a  tool  is  now  available  that  can  perform  arithmetic  very  rapidly. 

The  CRAY-1  computer  (Reference  6)  can  perform  over  100  million 
additions  in  one  second.  Previous  aerodynami cists  had  no  such  tool  and 
their  technology  developed  along  other  routes.  Today,  we  discretize 
partial  differential  equations  into  finite  .difference  algebraic 
relationships  where  arithmetic  can  be  used.  Since  the  large-scale 
computer  is  very  good  at  this,  the  approach  is  successful. 

We  are  presently  developing  methods  to  exploit  the  computer 
capabilities.  This  new  field  is  called  "Computational  Aerodynamics . " 
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SECTION  II 

GOVERNING  EQUATIONS  IN  AERODYNAMICS 

Various  levels  of  sophistication  can  be  employed  to  attack  a  problem 
in  aerodynamics.  The  approach  one  takes  depends  upon  the  accuracy  required, 
time  and  funds  available,  etc.  An  attempt  has  been  made  to  catalog  the 
prediction  methods  to  help  standardize  the  nomenclature.  Figure  1  lists 
the  methods  and  restrictions  while  Figure  2  graphically  shows  the  limits 
of  the  methods . 

In  this  section  the  governing  equations  of  fluid  mechanics  (Navier- 
Stokes  equations)  will  be  derived.  Prior  to  this  derivation  however, 
some  vector  analysis  (Reference  7)  relationships  must  be  reviewed. 

1.  VECTOR  ANALYSIS 


The  following  definitions  will  be  required, 
coordinates  are  given. 


Examples  in  Cartesian 


Vector 

Elemental  Area 

Dot  Product 
Divergence  of  Vector 

Curl  of  Vector 


Dyadic 


v_=j_u  +j  v+kw 

dA=|dAx+jdAy  f  kd  Az 

V-dA=u  dAx  +  vd  Ay +wd  Az 

v.v=  <3u  +dr  +  dw 

~  dx+dy  +  dz 


VxV  = 


I  1  J< 

d_  d_  ± 

ax  ay  az 


P=i 


u  V 

w  i 

cr 

r 

r 

II 

12 

13 

r 

0" 

T 

21 

22 

23 

TV 

r 

cr 

31 

32 

33 
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Dot  Product  of  Dyadic 

V-P=i(u|oj  +  Vt2|  +  wt3|  ) 

+1(UTI2+V,r22+Wr32) 
+  *(ur|3+VT23+  w<r33> 


Divergence  of  Dyadic 


dr  dT 
r21  +  31  \ 

dy  dz  1 


,,/arl2+a°22+aT32\ 
1\  dx  dy  dz  / 

.dr.,  (3t-__  . 

+k  (-t'it1  ^r) 

ox  ay  dz.  / 


Substantial  Derivative 


DV  dV  dV  w2 

5ri7  +  <¥-v)y=  gf  +v1-Vx(vxv) 


<#>V-dA=///(V-y)dV 

Green's  Theorem 

2.  DERIVATION  OF  NAVIER-STOKES  EQUATIONS 

The  governing  equations  of  fluid  mechanics  (References  3  and  9)  are 
derived  from  statements  of  the  conservation  of  mass,  momentum  and  energy 
for  an  arbitrary  control  volume. 
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3.  CONTINUITY 

Statement  of  the  Conservation  of  Mass 

Net  Outflow  of  Mass  _  Decrease  of  Mass 

Through  Surface  in  Control  Volume 

But  Green's  Theorem  states 

#>/>V*dA=/j7(V;AV)dV 

Hence,  after  substitution 

///fir  +v-/»y]dvo 

Since  the  control  volume  (V)  is  completely  arbitrary,  the  integrand  of 
the  integral  must  vanish. 


dp 

at 


+V-/3  v=o 


Continuity  Equation 


4.  MOMENTUM  EQUATION 


Statement  of  the  Conservation  of  Momentum 


1 )  '  •  ' 

■Sum  of  External 

Rate  of  Change 

Forces 

of  Momentum 

F 

= 

rrr  DV 

fJIPD fdV 

F 

sum  of  stresses  x  Area  in  unit 

vector  directions 

= 

<#>P-dA 

e  P 

— 

stress  dyadic 

or  stress  tensor 
or  stress  matrix 
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A  dyadic  is  treated  as  a  "double  vector"  and  is  manipulated  as  such.  It 
is  symmetric  and  composed  of  three  normal  stresses  and  six  shear  stresses 

•'  Ps(-p+\v*y)l+^(vy+yv) 

Green's  theorem  is  now  used  as  follows:  . 

F  =  <^>p.dA=///(V.P)d¥ 

Hence, 

//;[/>5f-v-P]dv=o 

Similarly,  the  integrand  must  vanish.  Momentum  Equation 


5.  ENERGY  EQUATION 

Statement  of  the  Conservation  of  Energy 


Rate  of  Heat 
Added 

+ 

Rate  of 
Work 

Done 

Rate  of  Change 
-  in  Internal  Energy 

dQ 

dt 

+ 

dW 

dt 

= 

dE 

dt 

<$p  dA 

¥ 

f  •  y 

3 

xcr^dv 

^[^H-P-vj-dA 

+ 

///*■§?  dY 

where 

e=  CvT  +^*  and  <^_=  kVT 

Green's  Theorem 

///V'[p-y+^]d¥  =//J>^dV 
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Similarly,  the  integrand  must  vanish. 

^§7=v-[|>.v  +  i] 


6. 


or 


Dt 

DIVERGENCE  FORM  OF  EQUATIONS 

Add  the  product  of  ’ e '  times  continuity  to  the  energy  equation 

e  (ff+v'^)+',Dt a7+^9t+8V'/,y+/,-'Ve 

=  +v-pv e  =V.jp-V  +  |] 


-^y.(pe)  +  V«[pVe- P-V- cj_]  s  0 


Similarly,  adding  the  product  of  times  continuity  to  the  momentum 
equation  produces  the  divergence  form  of  the  momentum  equation. 


^(/>v)+  v*[p  v  v-p]  =  o 


Note  that  the  conservation  of  mass,  momentum  and  energy  can  now  be  written 
in  identical  form  using  the  divergence  vector  in  Cartesian  coordinates 
(References  10  and  11). 

£LT+iI  +  if:+iiL=o 

dt  dx  dy  dz 


where 


U  = 


p 

su 

pa 

,u2-ctm 

pv 

i  E  = 

?u  v-<T|2 

pw 

?uw-r)3 

Pe 

Sue-ucTjj-  vt,2”wt]3~  kTx 

,etc. 
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wherein  -  _p-r\V-V+ 2/iux;  rj2  =  ft(uy+vx);r|2=At(  uz+wx) 


The  three  equations,  (2  scalar,  1  vector),  contain  four  unknowns,  i.e. 

(V_,  p,  p,  T).  The  equation  of  state  is  needed  to  close  the  system.  For 
the  present  case,  an  ideal  gas  is  assumed. 

p  s  PR  T;  equation  of  state. 

This  is  the  fourth  required  relationship.  Values  for  the  transport  and 
thermodynamic  properties  (u,  k,  Cv,  R)  are  also  required  as  input.  These 
equations  with  appropriate  boundary  conditions  are  capable  of  representing 
nearly  any  aerodynamic  problem  in  the  aviation  field. 
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LEVEL  TYPE 


LIMITATION 


COMPLEXITY 


COMPUTER  TIME 


0 

EMPIRICAL 

QUALITATIVE 

AL6EI 

BRAIC 

1 

LINEAR 

SMALLtfM^I 

II 

IINVISCID 

NO  SEPARATION 

} 

r 

III 

NAVIER -STOKES 

NO  RESTRICTION 

PARTIAL  DIFF  EON 

SEC 


t 

HOURS 


Figure  1.  Aerodynamic  Prediction  Methods 
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SECTION  II.I 

SURVEY  OF  AERODYNAMIC  PREDICTION  METHODS 

In  this  section  lower  forms  or  approximations  of  the  Navier-Stokes 
equations  will  be  derived. 

1.  Level  III.  Navier-Stokes  (References  12  through  20) 

iy  +  ii  +  if+iG  n 
at  ax  ay  dz  'u 

a.  Level  III. A  Parabolized  Navier-Stokes  (References  21  and  22) 

aus0  •  a_E  __  a_F  _  a_G 


at 


ax  ay  dz 


where 


E= 


^2 
l fiU  +P 

p  uv 
P  uw 
/jue 


Viscous  terms  only  in  F  and  G 


b.  Level  III.B  Two-D  Boundary  Layer  (References  3  and  23) 


P  = 


-P  /xlly 
/j.Uy  —  P 


11 


AFWAL- TR-82 - 3031 


(/jU^+f/av^O 

(AU2+P)x  +  (/JUV-/zuy)y*0 

Py  *  0 

(/ouH)x  +  (/>vH-ur-kTy)y=0 


2.  LEVEL  II 

a.  Level  II. A.  Inviscid  (Euler)  (Referebce  24) 


fjL* 0  )  k=0 

P  — pi 


V-PV  =  0 
v-(/>v  v+pl)=o 
V-[pVH)*Q 

p 

where  H  =  e  +  — 


By  combining  the  last  equation  with  the  first,  one  finds  that  H  =  constant. 
Hence  one  differential  equation  reduces  to  an  algebraic  equation  thereby 
reducing  the  computer  time  to  solve  the  system. 

Alternate  forms  of  the  momentum  and  energy  equations  are  often  used. 


DY  Vp  dV^v2 

57'~  T 1  J7  + VT 


-  Vxw  and 


,L iv2 
2 
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b.  Level  II. B.,  Inviscid,  Irrotational  (Full  Potential)  (Reference  8) 

One  of  the  greatest  simplification  arises  when  the  vorticity  in 

the  flow  field  is  zero,  to  =  0.  This  implies  that  the  viscosity  vanishes 

(which  was  already  assumed)  and  no  shocks  exist.  In  practice  this  means 

P  o 

that  Mn  <  1.5  (for  which  a  1%  total  pressure  drop  occurs),  and  (Mn  = 

1.5)  <  2.46  or  only  weak  shocks  are  permitted.  When  this  occurs  a  velocity 
potential  can  be  introduced. 

v=  v<£ 

Automatically  this  insures  that  the  vorticity  vanish. 

cu  =  VxV  =  VxV<£  =  0 

since  the  Curl  of  the  Gradient  vanishes  identically.  The  governing  equations 

becomes  as  follows:  _  _  ,  ^ 

V-pV4>  =  0 

?PV(V<p)Z  +  V{paZ)--  0 

and  a2=a2  -  (V<£)2 

By  eliminating  P  these  equations  produce  the  full  potential  equation. 

a2  V2<£  =V<£-(V<£  -VV<£) 

Expanding  this  equation  in  Cartesian  coordinates 

a2(<^xx  +^yy  ^xx  +  V^yy  +  W^zz 


where  U  =  <£x 

V  =  *y 


+  2UV  4>xy  +  2  U W  <pxz  +  2VW0yz 
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3.  LEVEL  I.  LINEARIZED  EQUATION  (REFERENCES  8  AND  25) 

Level  I  is  further  restricted  by  simplifying  the  non-linear  partial 
differential  potential  equation  to  make  it  linear  for  which  analytic 
methods  of  solution  have  been  developed. 

Assume  small  perturbations  (Reference  8) 

U*V00  uU'(x,y,2) 

V=  V'(x,y,z) 

W=  W/(x,ytz) 

or  <p=XTJoo-b^/ 

Hence-(l-M^)  <^'x  +  ^y  +  ^2«o 
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TABLE  1.  REQUIRED  NUMBER  OF  BOUNDARY  CONDITIONS 
Level  III.  Navier-Stokes 

Initial  Cond.  B.C. 


Variabl es 

t 

X 

,y 

z 

u 

2 

2 

2 

V 

1  , 

2 

2 

2 

w 

1 

2 

2 

2 

p 

1 

1 

1 

1 

T 

1 

2 

2 

2 

TOTAL  27  +  5  =  32 

Level  III. A.  Parabolized  Navier-Stokes 

Initial  Cond. 


Variabl es 

D 

y 

Z 

u 

1 

2 

2 

V 

1 

2 

2 

w 

1 

2 

2 

p 

1 

1 

1 

T 

1 

2 

2 

TOTAL  18  +  5  =  23 


LEVEL  III.B.  Two-Dimensional  Boundary  Layer 


Variabl es 

X 

y 

u 

1 

2 

V- 

0 

l 

P 

i 

i 

T 

l 

2 

TOTAL  6+3 
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TABLE  1.  REQUIRED  NUMBER  OF  BOUNDARY  CONDITIONS  (Continued) 

Level  rr.A.  I'nviscid  (Euler) 

Variabl  es 

u 
v 
w 

P 


'Level  II. B.  Inviscid,  Irrotational  (Full  Potential) 

Variables 

♦ 


Level  I.  Linearized  Equation 

Variables 

♦ 

TOTAL  =  6 


X 

y 

Z 

2 

2 

2 

TOTAL 


X 

y 

z 

m 

m 

■ 

1 1 

■  I 

■ 

i 

i 

■ 

■ 

TOTAL  =  12 
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TABLE  2.  SUMMARY  OF  BOUNDARY  CONDITIONS 


Level 


Title 


No.  Terms 


No.  B.C, 


Restriction 


Navier-Stokes 

A.  Par.  N.S. 

B.  Boundary  Layer 

Inviscid 


None 

2 

~  =  0 


—  «  1 
u 


A.  Euler 

B.  Potential 


u  =  0 

weak  shocks 


Li  near 


U  =  V  +  U' 


TABLE  3.  DIAGRAM  OF  AERODYNAMIC  PREDICTION  LEVELS 


P.N.S. 


*  «  1 
u 


u  =  0 


Invi sci d 
(Eul er) 


oj  =  0 


~  +0 

3t 


Full 

Potential 


Li  near 
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SECTION  IV 

ANALYTIC  SOLUTION  OF  BOUNOARY  LAYERS 

One  of  the  first  branches  of  fluid  mechanics  to  exploit  numerical 
methods  was  the  boundary  layer  field.  Most  of  the  important  features  of 
CFD  can  be  demonstrated  by  examining  boundary  layer  solution  techniques. 


Consider  a  steady,  incompressible,  2-0  boundary  layer  flow  (Reference 


26). 


ux  + 


0 


uux+vuy=uex+yUyy 

with  boundary  conditions 
u(x,o)=0 

v(x,o)=0  u(o,y)=U<~ 

u(x,oo)=  eu(x) 

For  simplicity,  consider  Blasius  flow  in  which  Ue  =  U  . 

co 


A  transformation  is  used  to  realign  the  coordinates  to  obtain  better 
numerical  accuracy  and  to  simplify  the  boundary  conditions  (Reference  27). 


Let 


V*«klfc+’»*F] 
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Using  the  chain  rule 


jL.  jL  +  t> 

drj  d(  ^ 


<3^7 


d  g  d.  .  d  d 

Ty~-*y7i+1yr,^  va7 


The  transformed  equations  become  as  follows: 

Vt?+  F=-2£f^ 

F__-  V  F  =  -  2.£  FR. 
VV  V  i 


The  right-hand  sides  of  these  equations  are  zero  for  the  situation  where 
F  =  F(n).  This  type  of  flow  is  called  "similar"  in  that  all  velocity 
profiles  at  different  stations  can  be  collapsed  onto  one  similar  curve. 
The  governing  equation  then  becomes  an  ordinary  differential  equation 
which  was  first  solved  by  Blasius  (1908)  using  an  infinite  series  (Refer¬ 
ence  26) . 

Let  V  =  -f.  Hence  the  governing  equation  for  Blasius  flow  becomes: 


f'"+  ff"=  0 

with  boundary  conditions 

f(0)=0 

f'(0)=0 


A  series  solution  is  obtained  as  follows: 

Assume  f  ■  :  °o  +  °,  v  +  ^  +  a3  ^  +  a4|?+  ... 

The  solution  obtained  by  substitution  of  this  series  into  the  governing 
equation  and  applying  inner  boundary  condition  is: 

a2V2  a27?5  a2H^?8 

2!  5!  +  8! 
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And 


f//(0)='0.46960=a2 


is  obtained  from  the  outer  boundary  condition. 


The  friction  coefficient  at  the  wall  is 


2  /*■  (du  \ 

^2  Uy  Jo* 


Z^1  F.  (0) 


u 


oo 


or 


„  _  0.66411 

'"“vS" 
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SECTION  V 

NUMERICAL  SOLUTION  OF  BOUNDARY  LAYERS 

With  the  advent  of  the  digital  computer,  numerical  techniques 
(Reference  3)  were  developed  because  they  removed  all  the  limitations 
required  to  obtain  analytic  solutions. 


The  same  transformed  boundary  layer  equations  and  boundary  conditions 
are  used  for  either  the  analytic  or  numerical  approach. 

V^-F  ;  V(0)=0 

F^-VF^O  -  F(0)=0  and  F (<*>)  =  I 


The  velocity  profile  is  discretized  into  a  series  of  points  at  equal 
intervals  in  n,  i.e.  An  =  constant  (Figure  3).  By  using  Taylor  series, 
relationships  between  neighboring  points  can  be  obtained  (References  28, 
29  and  30) . 


C  -  C  J.  C7  A  .  r-"  Al7  ,  C///  i  r-,V 

Fn+r  F„+FnA,'  +  F„  2l  +  Fn  3T  +  Fn 


A V 
n  4 ; 


Similarly 


AV2 


Fn _ ,  =  F  -  F'  At,  +  F p'"— -i-  +  F 

n  1  n  n  /  n  2!  n  31  r 


n  4 1 


By  substracting  these  two  relations,  one  can  obtain  Fn . 

F'  Fn+I-Fn-I  r/^A7?2  , 
n"  ZAV  "  6 

By  adding  the  two  relations  one  obtains  Fn . 


p'/_  ^n+ 1-^^n  +  Fn-i  _  p iv  At/ 
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2 

If  we  truncate  the  original  Taylor  series  after  the  An  terms  we  obtain 
an  algebraic  finite  difference  relationship  for  both  the  first  and  second 
derivatives.  Since  the  governing  boundary  layer  equations  have  only 
second  derivatives,  this  "second  order"  method  is  sufficient  for  our 
purposes . 


Using  these  finite  difference  relations  the  governing  equation  becomes 

~an  ^n+l +  bn  ~  cn  Fn - 1  3  ^ 

b"*2  A 
At7 

cns  *'+ vn  “g" 

where  is  obtained  by  numerical  integration. 

V=  -/Vd-rj 

o  ' 

Using  the  trapezoidal  rule 

Vn=-[F|  +  2F2+2F3+-2Fn-l  +  F„J 

or  Simpson1 s  Rule 

V -[F, +  4F2+ 2F3+  - 2Fn-2  +  4F„-,  +  t,]  ^ 


A  7} 

T 


Therefore,  at  each  point  in  the  field  the  governing  equation  can  be  ex¬ 
pressed  only  in  terms  of  information  available  at  adjacent  points. 
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The  equations  for  each  point  developed  a  regular  pattern. 


Written  in  Matrix  form,  this  is  called  a  tri -diagonal  matrix. 

[a]  IF!  =  IBI 


where 

A  = 

We  have  N  linear  algebraic  equations  with  constant  coefficients  (for  each 
iteration  cycle)  which  contain  N  unknown  values  of  F.  To  solve  this  system 
we  might  consider  using  the  standard  Cramer^s  Rule. 

F  =  IA-BI 
n  — F 
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This  computation  requires  (N+l)!  multiplications.  An  estimate  of  the 
required  computer  time  can  be  made  for  the  CRAY-1  which  accomplishes 
80  million  multiplications  per  second. 


N  (N+l)! 

3  24 

10  39, 916., 800 

18  1 .2  x  1017 


CRAY  Time 

3  x  IQ'7  sec 
0.5  sec 
47  years  ! ! ! 


The  surprising  escalation  of  computer  time  with  the  number  of  grid 
points  shows  the  impractical ity  of  using  Cramer ks  Rule.  Fortunately  a 
simple  algorithm  exists  for  solving  tri-diagonal  matrices  (which  is  a  form 
of  Gaussian  elimination)  and  is  commonly  known  as  the  Thomas  Algorithm 
(Reference  3). 


Thomas  Algorithm 


-a_  F_  ,,  +  b„  F„ -c_  F  |  =  d  where  a  >0 
n  n+i  n  n  n  n-l  n  n 


Assume  the  existence  of  a  linear  relationship 

=  en  ^n  +  l  +  fn 


»n>0 

cn>0 

bn-  an  +  cn 


Hence 


Fn-r  9n-t  Fn+fn-l 


Substituting  this  relationship  into  the  original  tri-diagonal  and  solving 
for  Fn  produces  the  following: 


F  =/fn _ \  f  4. .  +  /dn+  cn  fn-l \ 

n\bn-cnen_|/  n  Vbn-  c„  en_  ,/ 


n  n  n  - 1‘ 
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which  confirms  the  original  assumption  of  a  linear  relationship, 

F  =  e  F  ,  +  f  . 
n  n  n+1  n 


Therefore 


"  Wn-I 


Wn-I 
"bn-cn9n-l 


f.= 


The  solving  procedure  is  then  quite  simple.  n>l 

Starting  at  the  surface  where  F-|  =  0  then  ^  =  e.j  Fg  +  f-|  3  0 
We  see  that  e-j  =  0  and  f-j  =  0  (since  F^  is  arbitrary  in  general). 

At  the  next  point 


and 


Vb3’C3e2 


d3+C3f2 
3  ^3~c3®2 
etc. 
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This  sweep  procedure  is  continued  until  e^  -|  and  ^  are  obtained.  At 
this  point,  we  sweep  back  down  and  evaluate  F  ,  using  the  fact  that  F^  =  1  . 

p  s  e  F  +  f  *  e  +  f 
N-l  N-l  N  N-l  N-l  N-l 

and  using  previously  evaluated  en  and  fn> 


F  ~  &  F  +  f 
N-2  N-2  N-l  N-2 


.  etc. 

Continue  until  all  F  are  found. 

n 

This  procedure  is  efficient  and  accurate  in  that  round  off  errors 
are  seldom  encountered.  Only  three  multiplications  and  two  divisions  are 
required  at  each  point.  Hence,  the  computer  time  is  proportional  to  only 
kN,  which  far  surpasses  the  efficiency  of  Cramer's  rule. 


Figure  3.  Discretized  Velocity  Profile 
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1.  SUMMARY  OF  NUMERICAL  PROCEDURE  FOR  BOUNDARY  LAYER  CALCULATIONS 

By  observing  the  numerical  procedure  and  strategy  employed  in  the 
solution  of  boundary  layers,  we  may  learn  some  lessons  that  will  be  useful 
in  solving  the  Navier-Stokes  equations.  The  following  are  some  important 
steps  in  the  process: 

1.  Formulate  the  governing  partial  differential  equations.  Insure 

that 

Number  of  Unknowns  =  Number  of  Equations 

Number  of  B.C.  =  Order  of  Highest  Derivative 

2.  Transform  Governing  Equations.  This  simplifies  boundary  condi¬ 
tions  and  achieves  better  numerical  accuracy. 

3.  Convert  to  Finite  Differences 

4.  Employ  Proven  Solving  Scheme  compatibl e  with  computer. 

5.  Satisfy  Stability  Constraint 

6.  Keep  it  Simple  Stupid  (KISS) 

Maintain  lowest  order  of  derivatives  in  system  of  equations. 

Avoid  elegance  and  sophistication. 

Make  the  computer  do  the  work.  Minimize  your  work  off  the 
computer.  (This  is  frequently  opposite  to  the  strategy  in  analytic  efforts.) 

7.  Integrate  numerically,  not  analytically.  Use  graphics  to  minimize 
print-out  and  achieve  data  compression. 
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SECTION  VI 

TRUNCATION  ERROR  ANALYSIS 

Let's  return  to  the  formulation  of  the  boundary  layer  analysis  but 
include  higher  order  terms  (References  31,  32,  33  and  34) 

F-VF--gn  Fn_|+&„  Fn-cn  Fn+rdn  =  0 

At)2 

where  a  ,  b  ,  c  retain  the  same  previous  definitions  but 
n  n  n 

dn=  ( FIV-  2VF"')  A-ry4  =  <j>  A^4 
12  12 

This  term  is  representative  of  the  truncation  error  in  approximating  the 
differential  equation  by  finite  differences. 


An  estimate  of  the  error  in  wall  friction  can  be  obtained  by  integrating 
the  governing  equation. 


CO  CO  GO 

/r/dn-/VF'dn=/ 
0  0  0 


t 

’.fn-i+Vir 

~yor 


cn  Fn+I 


drj 


CO 


/<£<*  V 


CO  CO  2 

f<pd7j  -2f  [?')  6-rj  .7 
o  o 
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Hence 


A£(0)  =  A7?2  ( .7 )  a  +  .  06  At?2 

F'(O)  12  (.4696) 


A  simple  relationship  for  the  error  in  wall  shear  stress  has  been 
obtained  by  including  the  next  term  in  the  Taylor  series. 

1.  RICHARDSON'S  EXTRAPOLATION 

Numerical  experiments  can  be  conducted  in  place  of  evaluating  the 

function  o/°°  <}>dn  for  more  complex  problems.  Recognizing  that  the 

truncation  error  term,  for  second  order  methods,  will  be  proportional  to 
2 

An  ,  a  technique  can  be  developed  which  can  extrapolate  the  results  to 
zero  error. 

Write  two  experssions  for  the  error  for  two  different  step  size 
calculations 

Fj  -  F^exact)3  k  At 
F  (exact)  =  k  A772 


By  dividing  these  two  equations,  k  (which  in  general,  is  difficult 
to  evaluate  analytical ly )  can  be  eliminated.  The  exact  value  then  can  be 
predicted. 


/  Ai7i  \2 
FW.)  =  !•'-(— ) 


/ 


This  is  known  as  Richardson's  extrapolation  and  is  a  practical  method  for 
error  estimation. 
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Below  is  a  tabulation  for  Blasius  boundary  layer  calculations  of  the 
error  term  in  wall  friction  for  different  An  step  size. 


TABLE  4 

TABLE  OF  ERRORS  FOR  BLASIUS  BOUNDARY  LAYER  CALCULATIONS 


Art 

F'(o) 

Error 

%  . 

HH9 

F'(o) 

Richardson '  s 
Extrapol ation 

Error 

0 

.46960 

0 

— 

— 

.1 

.46966 

+  .00006 

— 

.146 

.46959 

-.0001 

.2 

.46988 

.00028 

— 

.179 

.25 

.47004 

.00044 

.09% 

.180 

.46947 

-.00001 

.5 

.47140 

.00180 

.38 

.184 

.46959 

-.00013 

1  .0 

.47718 

.00758  ' 

1  .6 

.194 

.46957 

- . 00001 

2.0 

.50 

.0304 

6.5 

.194 

.330476 

.-.1391 

5.0 

1  .39 

.92 

200 

.94 
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SECTION  VII 
STABILITY  ANALYSIS 

-The  major  issue  in  choosing  a  finite  difference  algorithm  is  its 
stability  characteristics.  Of  the  many,  many  ways  to  represent  a  partial 
differential  equation  by  finite  differences  only  a  few  are  stable.  (Many 
come  forth  but  few  are  chosen).  The  best  way  to  appreciate  this  is  to 
consider  the  roots  of  a  very  high  degree  polynomial  of  order  n.  To  insist 
that  no  positive  real  root  exists  when  n  equals  several  thousand  is  a 
severe  restriction.  This  is  equivalent  to  demanding  that  no  positive  real 
eigenvalues  be  permitted  in  the  large  solution  matrix  representing  the 
entire  computational  grid. 

Since  many  points  exist  then  this  is  obviously  a  very  severe  constraint. 
Therefore,  most  of  us  should  use  only  proven  algorithms  and  not  invent 
new  ones . 

1 .  MODEL  EQUATIONS 

To  study  stability,  we  shall  use  linear  model  equations  to  demonstrate 
the  key  features  (Reference  35,  36  and  37)  . 

The  boundary  layer  equation  has  terms  that  represent  advection 
(convection)  and  diffusion  of  the  following  form: 

u,+uux+vuy=„uyy 

These  two  processes  can  be  represented  by  two  simple  model  equations. 

Advection  U.  +  cU  =0  Simple  wave  equation. 

t  X 

Diffusion:  U  =  vU  Heat  equation. 

L  XX 

By  examining  the  stability  of  these  two  equations,  the  major 
limitations  can  be  assessed.  This  greatly  reduces  the  amount  of  labor 
i nvol ved . 
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2.  VON  NEUMANN'  STABILITY  ANALYSIS 

Let  us  first  consider  two  possible  finite  difference  representation 
of  the  simple  wave  equation  (Reference  31). 


1.  Forward  Time  -  Central  Space  (FTCS)  Differencing  Scheme 


U,+  CUX* 


n+l  n 
JJi  -Uj 


At 


c<Uj+|  ~  U]-  |)s  0 

2Ax 


2.  Backward  Time  -  Central  Space  (BTCS) 

t. 

t  =  nAt 

u.  feu  =— r- -re  - — 

f  x  At  2Ax 

C  Af 

Defining  a  ~  -r+  =  Courant  Number 
Ax 

Hence 


x  =  i  Ax  — ► 


t  i  i  i  i 


n+l  n 
Uj.Uj 


n  +  l 

ui+r 


n+| 

ui-i 


n+l  n  n  n 
Uj  — f  uI+l  +  u,  +  f  UUI 


2  wj+l T  UJ 


2  WJ- 


BTCS : 


n  n+l  n+l  n+l 

ur!ui+i+ui  -fui-i 
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Von  Neumann  (1944)  devised  a  method  for  analyzing  the  stability  of  finite 
difference  relationships-.  A  Fourier  series  expansion  of  the  solution  is 
performed  and  the  decay  or  amplification  of  each  mode  is  examined  to 
determine  the  stability  characteristics. 


Let 


UP=Un(t)e 


iax 


An 

where  U  (t)  is  the  amplitude  function  at  time-level  n  and  a  is  the  wave 
number .  (i=VH) 

Hence :  FTCS :  u"+  'eiax  =  ^  [- |eia  (x+Ax)  +  eiax  +|  eia  (x~Ax )] 

An+I 

and  Gain  =  G=  -  =  |  +  f  (e"iaAx-eiaAx) 

U  2 

or  G=  l-iasin  (aAx)=AMPLlFICAT10N  FACTOR 


For  stability,  | G j <  1,  since  the  solution  must  remain  bounded. 

-  2  2 
G  G=  1  +  c-  sin  aAx  <  I 

The  stability  criteria  indicates  that  cr  =  0  is  required,  which  is  im¬ 
possible  to  achieve. 

Therefore  FTCS  =  unconditionally  unstable. 

Likewise, 

un=un+l[l+!  (eiaAx-e"iaAx)] 

An+I  _l 

G=  —  =  [l  +  i<x  sinaAxj 


BTCS 

or 


GG= 


2  2 

l  +  cr  sin  aAx 


<1 
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IGI= 


-y/l+a-^sin^aAx 


For  stability  a  =  any  value.  Therefore  BTCS  is  unconditionally  stabl e 


3.  SUMMARY 


Simple  wave  equation 


FTCS:  any  a  or  any  step  size  is  unstabl e , 
BTCS:  any  a  or  any  step  size  is  stabl e. 

4.  EIGENVALUE  INTERPRETATION  OF  GAIN 

An  vt 
Let  U  =  Ae 


Hence  0"+,./WX#+A» 


An+|  \ a+ 

Therefore  G=JJ  =  =  1  +  XAt . . . 


0° 

G<  I  means  XAt  <0 
Since  At  >0  this  means 

\<0  or  no  positive, real  eigenvalue  in  time  is  permitted. 


5.  DIFFUSION  EQUATION 

Similarly  for  the  Diffusion  Equation 

vAt 


U  =  v  U  Let  d  =  — s  =  Diffusion  Number 


t  xx 


Ax' 


FTCS 


,  ,n  + 1  .  ,n  ,n  ,n  n  , 

U- -  Uj  v(Uj+|-2Uj  +Uj_|) 


At 


Ax' 
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Hence  the  Von  Newmann.  Stabil i ty  analysis  produces 

G=  I  -  2d  (I-  cosaAx) 


Si  nee 


Since IGI<  I 

1 2d  1(1 -cosaAx)-  I  j  <  I-,  But  (I  -  cos  a  Ax)  =  2  Maxvalue 
Therefore  d<-^--  Conditionally  Stable. 


Li kewise 
BTCS : 


G  =  l  +  2d  (I-  cos 


a  Ax)] 


IG  I  <  I  u>  Q. 

At>0  Always  true 

Ax2  >0 

d>0  Unconditionally  Stable 
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SECTION  VIII 
NUMERICAL  ALGORITHMS 

In  this  section  some  of  the  classical  differencing  schemes  will  be 
examined.  A  few  will  be  selected  which  are  representative  of  the  methods 
in  use  today.  Hopefully,  an  understanding  will  be  developed  from  studying 
these  few  cases  which  will  enable  the  reader  to  interpret  other  methods 
with  ease. 

The  cases  to  be  investigated  are  (Reference  31): 

1.  Leapfrog  Explicit 

2.  MacCormack  Explicit 

3.  Fully  Implicit 

4.  Crank-Nichol son  Implicit 

Roache  (CFD)  lists  several  other  methods  which  should  be  reviewed. 

1 .  EXPLICIT  METHODS 

An  explicit  method  is  one  in  which  all  of  the  values  on  the  right-hand 
side  of  the  difference  equation  needed  to  calculate  the  advance  time  level 
values  of  Un+^  are  known.  Methods  wherein  Un+^  also  appears  on  the  right- 
hand  side  are  called  implicit  and  generally  require  a  matrix  inversion  to 
calculate  a  new  time  level. 

2.  MIDPOINT  LEAPFROG  METHOD 

The  leapfrog  method  is  a  single  step,  second  order  accurate  explicit 
method  (Reference  36).  Since  the  method  is  central  in  time,  three  time 
levels  are  involved.  The  term  "leapfrog"  is  derived  from  the  fact  that 
the  new  values  are  calculated  at  each  other  time  level,  skipping  the  time 
level  in  between. 


36 


AFWAL-TR-82-3031 


Consider  the  wave  equation: 


du  c  <3u  _Q 
at  dx  ' 


n+l  n- 
U.  -  U- 

i  i 

2  At 

Recombi ni ng 


+  c 


Ul+|-Ul 

2At 


>r=  ur 


=  0 


-  <T  (  U 


CTCS 

o 


n  -U°  ) 
i  +  l  Ui-|J 


where  <x  =  cAt  Courant  number  (Ref  38) 
Ax 

The  stability  may  be  assessed  using  the  Von  Neumann  method. 

n  ad  ,  iax 
Let  U  =  U  (t)  e  ,  then 

fin  +  L  r,n“L  Ax  -iaAx^ 


U  =  U  —  cr  U  (e 
a  n+l  ^.n-l 


orlj "  'SU"  '  —  i(2o-sinaAx)On  =  aOn  +  UR 

Since  this  is  a  multi -time  level  method,  an  identity  relation  must  be 
added  to  determine  the  amplification  factor,  i.e. 


un=(i)un+o(un'') 


Hence 


An+I 

A,n 

U 

=  G 

U 

where  G  = 

a  1 

A,  n 

a  n-l 

1  0 

u 

U 

For  the  previous  one-level  method,  G  was  simply  a  number,  and  the 
stability  criterion  was  | G | <  1. 


For  the  present  case  where  G  is  a  matrix  G  -  X I  =  0 
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then,  X  <  1  is  the  stability  criterion.  Another  frequently  used  state¬ 
ment  is  that  the  "spectral  radius"  of  G  must  be  less  than  unity. 

a-\  |  " 

=  0 

I  (0-A) 
or  \2-a\-l=0 

Hence 

A=l(a±v/a2+4) 

Inserting  a  into  the  equation,  Roache  shows  that  X  =  1  for  a  <  1.  This 
i ndicates'that  the  leapfrog  scheme  is  marginal ly-conditionally  stable. 

3.  MACCORMACK  EXPLICIT  METHOD  (REFERENCE  4) 

An  extremely  popular  method  for  solving  compressible  flows  has  been 
the  method  of  MacCormack  (1971).  It  is  a  two  step  method  which  alternately 
uses  forward  and  backward  differences  (References  39,  40).  Although  each 
step  is  first  order,  the  result  after  the  two-step cycl e  is  second  order 
accurate  in  both  time  and  space. 


G-Xj  = 


Consider  the  model  equation:  1).  +  C  U  3  0 

t  X 


Forward 

Predictor  (.FTFS) 
Step 


u,n+I  -  U,"  +  c 
At 


Ax 


Backward 

Corrector  (FTBS) 
Step 


u!’+l-uin  +  cU^+2-U^2 


At 


Ax 


+cUi+rur° 

2Ax 


or  U 


n+l 

1 


+  U 


n+t 


cr(U,n+2 
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Performing  a  stability  analysis 

« ^  a. n  a n t 's  r.  .  .  *"  ioAxl 

2U  =  U  +  U  2  j(|-<r)  +  <xe  iaaxJ 


Ant,  An  iaAx 

andU  2=y  (|+cr— ere  ) 


.  .  An+i 

eliminating  U  c 


^n+l 


G®  “ : —  =  i  +4  (  l-cr  +<re  ia^x)(  l+cr-cre1^*) 
An  2  2 


or  G=  I- cr^(l-cosaAx)-i<rsinax 

Stability  Condition 

GG  =  l-(l-cosaAx)2<x2  (l-<r2  )<l 


This  condition  is  satisfied  provided  a  <  1.  Therefore,  MacCormack  ls 
method  is  conditionally  stable. 


4.  FULLY  IMPLICIT  (REFERENCE  31) 

The  methods  previously  described  are  explicit,  in  that  only  known 
values  at  previous  time  levels  are  needed  to  advance  the  calculation  to 
the  new  time  level  ^n+^.  We  will  now  discuss  implicit  methods,  which 
use  new  values  in  the  spatial  derivatives,  thereby  requiring  the 
simultaneous  solution  of  equations  at  (n+1)  in  order  to  advance  the 
calculations. 


Write  the  model  equation 

ut+cux=  0 

in  finite  difference  form 

using  FTCS  but  evaluating  the  advection  term  at  the  new  time  level  (n+1). 
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This  if  the  fully  implicit  method. 


un+l-u" 

J _ i_ 

At 


+c 


(un+l-Un+l) 

Vi+I  i-l,  / 


2Ax 


=0 


Error  (At,  Ax^) 


n-H  n_?(  n+l} 

ui  i  2vi  +  l  ui-l 


Using  the  Von  Neumann  stability  analysis 


aitH 

U 


(l+i<rsinaAx)=U 


Hence  G= 


1  +  io-sirraAx 


or  GG= - - <  | 

l+o-^sin^aAx 

Since  any  value  of  a  will  achieve  the  stability  condition,  the  fully 
implicit  method  is  unconditionally  stable. 


5.  CRANK  NICOLSON  IMPLICIT  (REFERENCE  41  AND  42) 

A  modification  of  the  above  implicit  method  is  to  use  FTCS  but 
evaluate  the  advection  term  at  the  average  between  the  (n)  and  (n+1) 
terms . 

For  the  model  equation  this  scheme,  developed  by  Crank-Nichol son 
(1947),  is  the  following 


un+l-un  (un+,-un+l)  (Un  -u"  1 

Ui  Ui- 1  tu1+l  ul-l)  „  „  ,  A  2  *  Is 

+  c  - - - +c - — - =  0  Error=  (At  ,  Ax^) 


At 


4Ax 


4Ax 
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The  amplification  factor  is 


aP+I 

U 


I  -  io^sinaAx 
2 

1  +  i^sinaAx 
2 


For  which  GG  =  1  which  is  marginally  stable. 
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-SECTION  IX 
SHOCK  WAVE  STRUCTURE 

Consider  a  traveling  shock  wave  in  a  long  tube  (Reference  8).  In 
most  aerodynamic  calculations  a  shock  wave  is  treated  as  a  discontinuity 
and  Rankine-Hugoniot  relations  are  used.  However,  a  shock  wave  in  nature 
has  a  continuous  structure  which  establishes  a  rapid  transition  from  one 
state  to  another.  To  analyze  this  situation  we  shall  utilize  the  unsteady 
Navier-Stokes  equations  in  one  spatial  dimension. 

Vt+Ex*° 


P 

Au 

whrere  V  = 

P  u 

;  E  = 

puz-% 

/3  u  e  -  ucrj  |  -  kTx 

dnd  e=  CvT~~  =  H— 

oj|--P+\  ux  +  2/xu x=-P  +  i /x  ux 

3 

2 

since  as-— ^  ancj  p,  Sfx(j) 


The  energy  equation  can  be  simplified  for  the  case  in  which  P^  =  =  3/4 


In  reality  Prandtl  number  for  air  is  0.72,  making  this  approximation 
quite  reasonable.  Using  this  condition  the  energy  equation  becomes 

P  4 

^ue-uo-j!  —  kTx=/ou(e+-^)--^yauux-kTx 


=/ouH--^/x  Hx  —  k  Tx  ( I  -y$-  Pr) 


First  an  analytic  solution  shall  be  obtained  and  then  the  numerical 
procedure  discussed. 
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1.  ANALYTIC  SOLUTION  (REFERENCE  43) 

First  transform  to  eliminate  the  variation  of  y(T)  effects. 

Let  d£=-p-  dx 


P 

pu 

4/a 

Hence  -j- 

pu 

+ 

p\£  +  P-u^ 

p,e 

t 

p u H-Hf 

=  0 
C 


For  steady  state  this  equation  may  be  immediately  integrated. 


/2U=  Cj 

pup-  +  R-u^.  =  Cj  C^=  C(  u+  P-u^ 
pu  H- Hf  =  C|C3=  C(  H-H^ 


The  last  equation  can  be  integrated  again. 


H  =  C  +  C  e 
3  4 


C 


Since  for  no  heat  addition  H  (5  ->•  ®)  must  be  bounded  we  conclude  that 
C^  ^  0,  and  hence  H  =  C^  =  constant  =  H  .  We  will  find  it  useful  to 
express  in  terms  of  the  acoustic  speed. 


y+l 

HT  2(y-l) 


thus 

_g_2_yp  _  x+|  2  r-|  2 

yu "  p  ~  2 ‘ °*  2  U 


Therefore  p  can  be  eliminated  from  the  momentum  equation. 
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2  2  2 
since  a  =  a  (u)  this  equation  can  now  be  integrated.  Inserting  a  into 

the  momentum  equation  and  regrouping  produces  the  following: 


n  u  . .  u  .  a- 
^  U.  ^  u.  "  yu 


uu 


NoterCj*^  Uj 


Vl r  ui <l+0) 

a=  2 
uf 


By  Integrating  once  more  the  distribution  of  u  through  a  shock  wave  is 
obtained.  _  , 

r+ 1  rz  p iuidx 


UwU  ,-a_^  d~a)  Ty  /) 


(l-~)(^-ara=C5e 


2.7  J  4/i 


The  final  integration  constant  Cc  is  determined  by  arbitrarily  setting 

b 

x  =  0  at  some  reference  point.  One  possibility  is  to  select  x  =  0  at 
u  =  a*  since  it  must  always  occur  in  the  interval  of  the  shock. 


Hence 


c5  = 


l-v'a" 


Mote  that  the  Rankine-Hugoniot  condition  is  recovered  in  the  asymptotic 
1 imit. 


or  U|U_  =  a*  (Prandtl's  Relation 

Reference  44) 
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2.  SUMMARY  OF  FIVE  BOUNDARY  CONDITIONS 

Given 


i  Hx=  0  atao 

Cg=  a  u  (0)  =  a*  t  u(0)*a* 

°  (va-a) 

3.  Numerical  Solutions  Returning  to  the  Governing  Equation 

V,+  Ex  =  0 

we  wish  to  conduct  a  numerical  integration. 

However  in  order  to  reduce  the  amount  of  programming  we  shall  make 
use  of  some  of  the  analytic  results.  Let  pu  =  C-j  and  H  =  and  hence 
only  the  momentum  equation  requires  numerical  integration. 


C|  =Vl 

C  .Vl+5  „+  \Z±i 

C2  p |  u,  "  (l  a)  2 7 


C  -  H  -Z+L 
C3“  Hl  “  2(t-I)  * 

C,=  Oor  H  =H. 


(^u)t  +  (pu2-  0jj)x  =  0 
u  [/°t+  (/>u)J  +/3ut  +/3UUX-  {<T|,)X  =  0 


Divide  by p 

I  Per  4 

ut  +  uux  +  ^(— “T-“UX>X*° 


Rewriting : 

uf  +  uGx*° 
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Where 


G 


=  u  + 


G(u) 


To  solve  the  above  equation  for  u(x)  we  shall  use  MacCormackls  two-step 
di fference  scheme . 


Forward 

Predictor 


Backward 

Corrector 


s  -L  u. 
2  i 


Care  must  be  exercised  in  evaluating  the  derivatives  in  G. 


8ACKARD  IN 
PREDICTOR 


FORWARD  IN 
CORRECTOR 


4.  FIVE  BOUNDARY  CONDITIONS 


\-P\  up  C| 

2and3.HI=C3=H2=1pfi54 

4.  u  {— x  )  s  u 

V  t  i 

5.  u(0)=a  Arbitrary  Origin  Reference 

or  ^2 
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•  SECTION  X 
ARTIFICIAL  VISCOSITY 


Before  obtaining  a  numerical  solution  one  should  examine  the  length 
scales  appearing  in  the  Navier-Stokes  equations.  To  accomplish  this  we 
shall  examine  the  dimensionless  independent  variables.  They  are  of  the 
following  form: 


We  first  see  a  time  scale  of  jj.  This  term  is  defined  as  a  character!' Stic 
time  (t  which  equals  the  time  required  for  a  particle  to  traverse  the 
computational  domain  (L).  A  characteristic  time  (t  is  a  good  measure 
of  the  time  for  transient  phenomenon  to  occur.  Generally  the  inviscid 
field  requires  about  3  t  to  attain  steady  state  based  upon  both  shock 
tunnel  and  numerical  experience. 


In  space  we  have  a  scale  length  of  L  which  is  derived  from  the  boundary 
conditions  imposed  at  the  edge  of  the  computational  domain.  The  other 
length  scale  is  which  is  proportional  to  the  mean  free  path. 

\=  I.62j/=  mean  free  path  —  10  ®  ft  at  sea  level 
a 

Since  in  practical  problems  these  two-scale  lengths  are  orders  of  magnitude 
apart  it  is  apparent  that  numerical  difficulties  should  be  anticipated. 


A  derived  intermediate  scale  ari-sing  in  solving  these  problems 
results  from  a  combination  of  the  previous  two  scales. 


S=yGT=  ==  =  Boundary  Layer  Thickness 
u.L 
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In  engineering  practice  we  attempt  to  honor  both  L  and  <5  but  disregard  x 
as  being  unimportant. 

Assume 


The  equations  and  numerical  solving  technique  are  unaware  of  our 
intentions  and  regard  the  small  scale  lengths  as  introducing  mathematical 
"stiffness"  into  the  problem.  Numerical  instabilities  occur  if  calcu¬ 
lations  are  attempted  with  this  disparity  existing.  In  computing  shock 
waves,  for  example,  oscillations  occur  (Gibbs  phenomena)  since  the  shock 
thickness  is  far  less  than  the  step  size  used  in  engineering  practice. 

To  eliminate  this  problem  an  aritifice  is  required.  Since  we  cannot  make 
Ax<X  and  solve  any  practical  problem,  let's  multiply  X  by  a  factor  (s)  to 
make  it  as  large  as  Ax. 

Let  Ax=/3X 

Therefore /3X=  1.6/3^  =  Ax 

~~ap~ 

Let  \i'  =  Bu  which  is  artificial  viscosity  added  to  the  equation  to 
remove  the  mathematical  difficulty  (References  45,  46  and  47). 


This  procedure  obviously  changes  the  physics  of  the  problem,  however, 
and  we  must  exercise  care  that  the  additional  viscosity  effects  are  no 
greater  then  the  truncation  error  in  the  finite  difference  scheme.  If 
this  objective  can  be  obtained  we  have  a  practical  solution  to  the  numerical 
difficulty.  Several  forms  of  artificial  viscosity  shall  now  be  discussed. 
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1.  NORMAL  STRESS  DAMPING  (REFERENCE  2) 

In  the  normal  stresses  two  viscosity  terms  appear, 

ojl  =  -  P+  X  V-V+  2ft  ux 
2 

where  X=- 

Rewrite  this  term  X=  +  -^-/3/x 

u 

whereas  x  =  Cell  Reynold's  Number 


This  scheme  has  been  used  successfully  by  McRae  in  treating  the  shock 
wave  for  a  cone  at  angl e-of-attack 

The  shear  terms,  e.g.  =  y  (Uy  +  vx),  are  unaltered  in  this  approach. 
Since  A  is  only  of  any  consequence  in  the  normal  stresses,  it  improves  the 
shock  capturing  capability  without  affecting  the  shear  terms. 

2.  VON  NEUMANN  RICHTMYER  DAMPING 

The  first  to  use  artificial  viscosity  merely  added  a  term  to  the  Euler 
equation  in  place  of  the  non-existent  Navier-Stokes  stress  terms,  e.g. 

2 

pu  +  P-p.  ux 
where  ft/  =  /oAx^ 

This  term  is  similar  to  a  turbulent  Reynolds  stress  and  is  of  second  order 
accuracy.  It  also  possesses  the  correct  sign  to  add  dissipation  to  the 
system. 
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3.  MACCORMACK'S  PRESSURE  DAMPING  (REFERENCE  4) 


MacCormack  rationalized  that  boundary  layers  possess  zero  normal 

pressure  gradient  and  therefore  an  artificial  viscosity  term  proportional 

to  ^  will  only  affect  shock  waves.  His  additional  damping  term  is 
3n2 

as  follows: 


[«“ii 


4.  UPWIND  DIFFERENCING  (REFERENCE  31) 

Some  differencing  schemes  possess  truncation  error  terms  that  behave 
like  artificial  viscosity.  The  upwind  differencing  method  has  this  feature 

Consider  the  following  model  equation: 


Construct  a  difference  scheme  that  is  central  in  time,  central  space,  for 
the  diffusion  term  but  upwind  for  the  convective  term. 


n+l  n-l  .  n  n  ,  n  A  n  n 

Uj  -Uj  (ui-Uj-i)  (ui+,-2Uj  +Uj_|) 

- TKT-  +  T* - *•' - 71 - 

Axc 


u  >0 


+  u- 


.  n  n, 

(ui+rui 1 


Ax 


u  <  0 


Meteologists  used  this  method  and  derived  the  title  "upwind"  for  the 
direction  bias  for  the  one-sided  differences. 


Let's  examine  the  truncation  error  based  upon  the  analysis  of  Hirt, 

nil  n 

Letuj  s  u.  ±  At  u-f +  -£- lift +’•■  for  time 

n  n  Ay2 

u.±|=u.  ±Axux  +  =2~  u^+’-for  space 
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Insert  these  relationships  into  the  difference  equation. 


Hence  IL+  UU  =i/U  ±U^U  »  +^Pwmd 

t  x  XX  2  XX  -Downwind 


orUt+UUx>[^  +  IUI^]uxx 


It  is  clear  than  an  effective  viscosity,  v  , 


Where  vA  -  v  I  + 


lUIAxl 


e  '  L'  Zv  J 


is  inadvertently  added  to  the  governing  differential  equation  by  the  finite 
difference  process. 


In  order  to  obtain  an  accurate  solution  it  is  clear  that  IJJj_Ax_< 


2. 
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SECTION  XI 

COORDINATE  TRANSFORMATION  PROCEDURE 

Very  soon  in  the  study  of  Computational  Aerodynamics  one  encounters 
configurations  which  cannot  be  described  by  a  Cartesian  Coordinate  system. 
Anal  tic  orthogonal'  coordinate  systems  exist  for  a  few  classic  cases,  i.e., 
cylindrical,  elliptical,  parabolic,  spherical,  conical,  paraboloid,  prolate 
spheroid,  oblate  spheroid,  etc.  However,  even  these  cases  are  certainly 
limited  in  application.  A  more  general  approach  is  required  to  analyze 
aircraft  components.  For  example,  consider  an  airfoil  of  arbitrary  shape 
(Reference  48) . 

Two  possible  grid  systems  may  be  considered:  (1)  Use  a  Cartesian 
grid  and  establish  an  interpolation  scheme  near  the  surface  to  describe 
the  boundary  condition  (Figure  4)  or  (2)  Generate  a  body-oriented 
coordinate  system  and  transform  the  governing  equations  (Figure  5).  The 
former  approach  retains  the  original  simple  form  of  the  governing  equations 
but  over  complicates  the  boundary  conditions.  In  addition,  thin  viscous 
layers  require  clustering  of  grid  points  near  the  surface  resulting  in 
further  difficulties.  The  later  technique  maintains  simple  boundary  con¬ 
ditions  but  adds  more  terms  to  the  governing  equations.  Clustering  of 
points  near  the  surface  is  readily  achieved  with  a  body-oriented  system, 
however,  the  task  of  grid  generation  is  an  added  burden. 

All  factors  considered,  the  body-oriented  system  is  a  clear  winner. 
This  method  shall  now  be  discussed. 

1.  BODY-ORIENTED  COORDINATES 

Consider  an  arbitrary  body-oriented  coordinate  system  (Reference  13) 

£  =  £(x,y) 

173'7(x,y) 

where  17= 0  descri bes  the  body  surface. 
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The  governing  equation,  Ex  +  Fy  =  0,  must  be  transformed  into  the 
new  coordinates.  To  accomplish  this,  the  chain  rule  of  differentiation 
is  used. 

Aaf  ±.  ± 

<3x  S  \  dr) 

±.c  _!+__! 

dy  y  d£  7  dij 

The  Jacobian  of  the  transformation  is 

^x 


The  transformed  governing  equation  becomes  (?  E  +  £  F  )  +  (nv  Eri  + 

x  s  y  5  x 

ny  Fn)  =  0. 

It  is  now  apparent  that  the  actual  functional  form  U,  n)  of  the 
transformation  is  not  required  because  only  the  derivative  metrics  (5  . 

X 

n  )  appear  in  the  governing  equation.  This  feature  facilitates  the  use- 
of  a  numerical  transformation  in  which  the  metrics  are  computed  by  finite 
di ffferences .  One  addition  operation  is  required,  however,  in  order  to 
maintain  a  simple  procedure. 


Let's  generate  the  metrics  by  means  of  the  inverse  transformation . 

*  =  *(£,  77) 

y=  y  (C,  77 ) 


£x=  Jy  V  Jy£ 

f  =  ~  J  X  •  77  5  JXf 

y  v  v  C 


and 


xv  yi 
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Using  the  inverse  transformation  metrics  an  alternate  form  for  the 
transformed  governing  equation  becomes 


J(y  E.-x  F.)+J(-yf  E +X.F  )  =  0 
V  Z  V  V  $  Tl  £  7J 


Dividing  by  J  and  manipulating  the  derivatives  produces 


(y<7E 


o 


+  xc  )s  0 

Zi 


A  A 

or  Eff  F  =  0 
Z  "7 


whereE  =  y_E-x  F 
1  1 


The  transformed  equation  is  now  identical  in  form  to  the  untrans¬ 
formed  equation,  however,  additional  terms  have  been  added  to  the  flux 
vectors.  The  main  reason  for  using  the  inverse  transformation  is  to 
facilitate  the  use  of  numerical  derivatives. 


Recall  that 


dS 

dx 


y  =  const. 


and 


77=const. 
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The  transformed  coordinate  lines  located  in  physical  space  readily 
permit  the  numerical  evaluation  of  x  but  not  5  ,  since  lines  of  constant 
n  have  been  identified. 


Ax 

YA£ 


77=  const. 


Ay 

Va^ 


£=const. 


Hence,  the  inverse  transformation  metrics  are  numerically  evaluated 
from  the  predetermined  grid  and  inserted  into  the  governing  equations. 

The  only  restriction  on  the  transformati on  is  that  it  be  one-to-one 
(single-valued)  and  the  Jacobian,  not  vanish  in  the  computational  domain. 
The  transformation  need  not  be  orthogonal  and  it  is  not  necessary  to 
evaluate  the  functional  form  of  the  transformation  (since  only  the  metrics 
are  required).  Also,  one  need  not  transform  the  velocity  components  which 
further  simplifies  the  procedure. 

Another  advantage  of  this  transformation  concept  is  that  equal  step 
sizes  can  be  employed  in  the  transformed  space  which  permits  the  use  of 
simple  finite  difference  operators  in  the  numerical  procedure.  This  is 
not  possible  in  the  interpolation  method  originally  considered  as  a  candi¬ 
date  . 

2.  CLUSTERING  OF  GRID  POINTS 

The  use  of  a  transformation  permits  the  contraction  of  grid  points 
in  regions  of  high  gradients.  Consider  a  function  E(x)  which  has  large 
values  for  the  higher  derivatives. 

The  gradient  E  expressed  by  finite  difference  is 

X 

_  ^i  +  l  —  ^i-l  _  E  Ax^ 

x=  2Ax  xxx  6 
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where 


N  =  number  of  grid  points  in  domain  L. 

If  E  =  L  (X)  ,  then  the  maximum  error  in  E„  becomes  (-Exxx  -7—  = 

•j—  n  X  0  iTlaX 

n(n-1)  (n-2 ) 

—  6  In 

The  percentage  error  is  shown  below  for  various  n  and  for  N  =  5. 

%  Error  in  Ex  (N=5) 


n 

%  Error 

i 

0 

2 

0 

4 

-4% 

6 

-13.31 

11 

-60% 

16 

'-140% 

Large  errors  result  for  high  gradients,  therefore,  one  might  conclude 
that  more  than  five  grid  points  are  required  to  reduce  the  error  to  an 
acceptable  level.  However,  another  approach  would  be  to  stretch  the  grid 
in  order  to  achieve  smaller  gradients  in  the  transformed  plane  thereby 
reducing  the  size  of  the  truncation  error  term. 

Let£=£(x) 

E  -f  f  ?*CEi+i-Ei-i  a£21 

*  *  e« '  ZA(  6 

Choosing  a  stretching  factor  of  the  form 

M-lxT 

Lm 
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now  produces  a  maximum  error  in  E  as  follows: 

X 


(-£x  E  )  max  =  n  ^IL"1  ^0. 

x  6  m  m 


6N 


where  N= 

A£ 

The  percentage  error  for  m  =  4  and  N  =  5  is  now  within  reasonable  limits. 


n 

Transformed  %  Error 

i 

-0.87% 

2 

-0.5% 

4 

0 

6 

+0 . 1 7% 

8 

0 

11 

-0.87% 

16 

-4% 

3 .  SUMMARY 

To  expedite  the  numerical  solution  of  flow  fields  over  arbitrary 
configurations  the  equations  are  transformed  into  a  body-oriented  coor¬ 
dinate  system.  ■  This  transformation  is  accomplished  numerically  and  points 
are  clustered  in  regions  of  high  gradients  to  minimize  truncation  errors  . 
In  addition,  equal  step  size  is  used  in  the  transformed  plane  to  simplify 
the  finite  difference  operators.  Only  the  independent  variables  are 
transformed  while  the  dependent  variable  velocity  components  remain 
oriented  to  the  original  (Cartesian)  system.  (It  is  not  necessary  to 
transform  the  velocity  components  unless  one  desires  to  eliminate  terms 
through  an  order  of  magnitude  analysis).. 

Also,  the  transformati on  need  not  be  orthogonal.  The  resulting 
transformati on  adds  a  few  additional  terms  to  the  governing  equations  but 
greatly  improves  the  accuracy  of  the  method  by  optimum  grid  positioning. 
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The  burden  of  the  method  is  therefore  placed  upon  the  Grid  Generation 
procedure  which  will  be  discussed  next. 


Figure  4.  Cartesian  Coordinates  with  Interpolation  on  Boundary 


Figure  5.  Transformed  Body-Oriented  Coordinates 
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SECTION  XII 

PARABOLIZED  NAVI ER-STOKES 

The  complete  Navier-Stokes  equations  offer  the  potential  to  solve 
any  problem  in  fluid  dynamics.  However,  the  procedure  is  the  most  costly 
of  any  prediction  method.  There  is  a  more  efficient  solving  procedure 
entitled  Parabolized  Navier-Stokes  (PNS)  that  can  be  used  under  some  con¬ 
ditions  (Reference  22).  These  conditions  occur  for  supersonic  flow  with 
no  streamwise  separation  (although  transverse  or  cross-flow  separation  is 
permi ssi bl e) . 

Under  this  physical  situation  the  elliptic  terms  in  the  x-direction 
(U  )  can  be  neglected.  No  downstream  information  affects  this  portion 

X  X 

of  the  flow.  For  this  situation,  a  great  mathematical  simplification 
arises  and  permits  the  use  of  PNS. 

We  shall  now  explore  the  development  of  this  method. 

Begin  with  the  complete  steady  Navier-Stokes  equations. 

+  *y  +  =  ® 


where 


E  = 


pu 

2 

pn  -cr.. 


/3UV-r 


12 


/3UW-T. 


13 


/>ue-uoj|-vT|2-WT|3-q 
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All  the  elliptic  terms  in  the  x-direction  are  contained  in  the  E 
vector. 

rn=-p+x(©+''y+^)  +  2'i® 


rl3^(Uz+©> 

vK© 


By  neglecting  these  first  derivative  terms  in  x,  we  obtain  the  PNS 
equations.  However,  in  practice,  all  viscous  terms  are  dropped  in  the  E 
vector  to  simplify  the  solving  procedure.  This  additional  simplification 
does  not  greatly  limit  the  method  much  more  than  the  original  assumption 
of  neglecting  only  the  U  terms. 

XX 


Hence 


\pu 


E= 


plfi  +  P 

/3UV 

/>UW 

/JUH 


=  Inviscid 


With  this  formulation  it  is  possible  to  march  in  space  (x-direction) 
in  a  manner  similar  to  marching  in  time  that  was  previously  utilized  with 
the  time-dependent  Navier-Stokes  method.  PNS,  however,  requires  up  to 
two  order  of  magnitude  less  computer  time  to  solve,  thereby  justifying 
the  approximation. 
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We  shall  now  demonstrate  the  method  by  investigating  the  flow  over 
a  body  of  revolution  traveling  at  supersonic  speeds.  The  axi -symmetric 
PNS  equations  for  adiabatic  flow  follow: 

Ex+-L(rF)r=H 


where 


pu 

E= 

/3U2+P 

;a2=a2+rHu2 
°  2 

P  uv 

r=A-(ur+  V 

crv 

r22=-P  +  XV-V+2/xvr 

F= 

cruv-r 

i 

2 

crv-c 

22 

Tee=-p+xv'Y+2^r 

0 

V-v  =  Ux+-j-(rw)r 

H  = 

0 

~°96/ r 

» 

By  using  MacCormack's  method  and  marching  in  x,  a  new  value  of  the 
E  vector  at  the  next  station  (x  +  ax)  is  obtained.  Resolution  of  the  E 
vector  is  required  in  order  to  obtain  the  primitive  variables  needed  in 
the  F  and  H  vectors. 

This  operation  requires  further  discussion. 

E  =  E.  +  Ax[H-4-(rFil 

2  1  L  rdr  J| 


Let 


pu 

A 

E2  = 

/>u2+  p 

= 

B 

P  uv 

2 

C 
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Therefore  v  = 

Au  A 


The  three  remaining  relationships  with  three  unknowns  p,u,P 

A  =  /JU 
B=/xi2+P 


2  =  /P  +  Iliu2 

ao  p  ^  2 


can  be  combined  to  produce  a  quadratic  equation  in  any  of  the  variables. 
The  variable  selected  for  resolution  is  Mach  number,  M,  a  combination  of 
all  three. 


r2  b2_  (i  +  yM2)' 


«2  2 
A  ao 


M2(I  +  ^M2) 


°rij2  /32-2/±£v£2- 2(7+1) 

M  =  - - - 

2  O'2-  -^-/32) 

The  positive  root  is  supersonic  while  the  negative  sign  produces  a 
subsonic  root.  A  predicament  in  this  solving  scheme  arises  in  that  a 
criteria  is  necessary  to  select  the  correct  root.  However,  a  more  serious 
limitation  is  encountered  in  that  the  subsonic  root  is  unstable. 


Recall  that  the  original  assumption  for  PNS  was  supersonic,  unsepa¬ 
rated  flow;  therefore  only  the  positive  sign  on  the  radical  is  selected 


MM2'2'  for  y  =  1. 4 

0.4  (9.8  -p) 


where  4.8</3  <  9.3 


(sec  Figure  6; 
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Care  must  be  exercised  in  selecting  the  first  grid  point  in  the 
boundary  layer  to  insure  that  it  remains  supersonic.  (Note:  Alternate 
procedure  have  been  developed  to  extend  the  method  by  setting  ^  =  0  in 
the  boundary  layer  for  M  <  1  and  eliminating  the  quadratic  root). 


Once  the  Mach  number  is  ascertained  the  primitive  variables  can  be 
determi ned . 


P= 


B 


\+yM‘ 


p-Qi  i+^m2) 

°o 


u  = 


B/A 


(1+l/yM  ) 


These  values  are  then  updated  and  the  marching  procedure  continued 
until  the  final  station  is  attained. 
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SECTION  XIII 
AIR  PROPERTIES 

Since  we  will  be  working  with  air,  a  review  of  its  properties  is 
appropriate.  In  particular,  we  will  require  numerical  values  for  the 
thermodynamic  properties  and  transport  properties. 

1 .  THERMODYNAMIC  PROPERTIES 

Air  is  a  gas  mixture  composed  of  about  80%  nitrogen  and  20%  oxygen. 
Traces  of  argon,  COg  and  H^O  vapor  do  not  appreciably  affect  the  thermo¬ 
dynamic  properties.  A  gas  mixture  can  be  treated  as  a  pure  gas  provided 
the  properties  are  evaluated  in  accordance  with  the  species  molecular 
weight.  Shown  below  is  a  table  of  the  individual  properties  of  the  air 
components . 

TABLE  5 

INDIVIDUAL  PROPERTIES  OF  AIR  COMPONENTS 


El ement 

Mol ecul ar 
Weight 

Mi 

Concentration 

C. 

R.  Gas  Constant 

ft2/sec2'  °R 

Yi 

n2 

28.02 

.7809 

1774 

1  .404 

32.00 

.2095 

1552 

1  .401 

39.94  . 

.0093 

1244 

1  .66 

Using  these  values,  the  air  properties  can  be  determined  as  follows: 

ZM.C.R. 

o  _  i  i  i 
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For  air,  the  required  thermodynamic  properties  are  the  following: 

Air,  M  =  28.966,  R  -  1716  ft2/sec2oR,  '(  =  1.40 

where  e=CvT,h=CpT,  and  p  =  yoRT 

ond.cv.jaj.cp.-za  , 

since  Cp/Cv  =  Y  and  C  p  -  R  +  Cv. 

These  therodynamic  properties  are  the  needed  values  to  be  used  in 
solving  the  Navier-Stokes  equations  in  the  region  where  air  behaves  as  a 
perfect  gas.  An  appreciation  for  the  limits  of  a  perfect  gas  is  required. 

2.  REAL  GAS  EFFECTS 

As  the  temperature  of  a  gas  is  lowered,  the  phase  will  change  from 
gas  to  liquid.  Further  decrease  will  solidify  the  liquid.  The  temperature 
values  depend  upon  the  pressure  level.  As  the  temperature  of  a  diatomic 
gas  is  increased  above  standard  sea  level  conditions,  vibrational  degrees 
of  freedom  arise  decreasing  y  and  accordingly  increasing  Cv  and  Cp. 

Since  the  molecular  weight  does  not  change,  R  remains  constant.  This 
domain  is  entitled  thermally  perfect,  calorically  imperfect.  As  temper¬ 
ature  is  further  increased,  dissociation  of  the  diatomic  gas  into  a 
monatomic  gas  occurs.  Higher  temperatures  cause  ionization.  All  these 
later  features  produce  departures  from  the  perfect  gas  law.  The  following 
is  a  list  of  typical  temperatures  for  which  the  changes  occur . 

°2  N2 
-360°  F  -346°  F 

-297°F  -321  °F 

5 ,000°F  1 0 ,000°F 

These  temperature  values  are  well  beyond  the  normal  limits  encountered 
in  low- speed  flight. 


Sol idi fies 
Liquefies 
Di ssociates 
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However,  at  hypersonic  speeds  this  is  not  the  case.  This  can  be 
demonstrated  by  computing  the  total  temperature  for  Mach  number  10. 

=—-*1+^—  M?  =  l  +  .2(IO)2  =  21 

tqd  * 

In  a  wind  tunnel  a  typical  value  of  To  is  2100°R.  Therefore, 

7^=  I00°R  =  360°F 

Hence,  1 iquefication  can  be  encountered  in  hypersonic  wind  tunnels  and 
indeed  must  be  guarded  against  to  avoid  erroneous  data. 


At  Mach  10  flight  speed,  T®  =  400°R  produces  a  To  =  8400  °R,  which 
clearly  dissociates  0^  • 

To  further  assist  in  acquiring  an  appreciation  for  the  various  regions 
for  which  real  gas  effects  are  encountered  (Reference  49)  a  map  of  the 
flight  corridor  is  presented  (Figure  7).  Note  that  during  reentry  both 
0^  and  dissociation  are  encountered.  In  this  range  the  thermodynamic 
properties  are  not  constant  and  hence  must  be  replaced  by  functional 
rel ationshi ps . 

use 

p  =  /5(h,p)  in  place  of  p  =  P/RT 
e=e(h,p)  in  place  of  e=  CvT 
T=T(h,p)  in  place  of  T  =  h/Cp 


In  the  governing  equations,  h  and  p  are  computed  and  interpolation 
tables  used  to  obtain  p,  e  and  T.  In  addition,  the  q  term  must  include 
the  heat  of  dissociation  to  account  for  recombination  heating  on  the 
surface.  The  net  result,  however,  is  that  perfect  gas  heating  is  nearly 
equal  to  real  gas  heating  due  to  the  fact  that  Lewis  number  is  near  unity 
implying  that  the  heat  transfer  by  diffusion  is  almost  the  same  as  by 
conduction.  Real  gas  effects  merely  redistribute  the  modes  of  heat  transfer 
without  changing  the  total  amount. 
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3.  TRANSPORT  PROPERTIES 

Three  transport  properties  exist  which  account  for  the  transport  of 
mass,  momentum  and  energy  throughout  the  gas.  The  relationship  for  these 
three  modes  of  transport  are  as  follows: 

ac. 

Mass:  m.  =/jD D=  diffusion  coefficient  (Fick's  law) 
Momentum:  r  =  ps  viscosity  (Stokes  law) 

_  .  <3T 

Energy:  k  k=  conductivity  (Fourier  law) 


To  solve  problems  in  fluid  mechanics  numerical  values  of  these  three 
transport  properties  are  required.  However,  in  practice  one  value  is 
prescribed  (u)  while  the  combinations  of  others  is  given  (k  and  D). 

Certain  combinations  of  these  coefficients  occur  naturally  in  the 
governing  equations  and  have  been  given  labels.  For  example. 

2  2 

Q(dissipation)  ry  /xV  /L  /xCp  /y  \ 

Q  (conduction)  kT/L  kT/L  k  \CpT/ 


uCp  2 

where  Pr~ — ; — =  Prandtl  Number-,  M 
r  k 


.j-fvi'l 

r- 1  VcpT/ 


Mach  Number 


Q(diffusion)  m;  Ah  /jOAhCi/L  ^DCp  /CiAh\ 
Q(conduction)  kT/L  kT/L  k  \CpT/ 


where 

(Lewis  Number )=  L 


/jDCp 
e=  k 
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Hence,  three  transport  properties  may  be  entrered  into  the  governing 
equations  as  follows: 


/i.=  (2.27xlO-8) 


3/2 


T+I98.6°R 


lb  sec 
ft2 


Sutherland's  Law 


,  _/>DCp 

L«snr 


1.4 
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SECTION  XIV 
BOUNDARY  CONDITIONS 

The  importance  of  the  boundary  conditions  in  solving  partial  differential 
equation  can  hardly  be  overstated  (Reference  50).  This  is  apparent  if  you 
note  that  the  only  difference  in  formulation  between  any  two  vastly  different 
flow  problems  is  the  location  and  value  of  the  boundary  conditions  since 
the  same  Navier-Stokes  equations  govern  the  interior  regions  for  all  fluid 
fl ow  probl ems . 

In  this  section,  various  boundary  conditions  will  be  explored.  First, 
the  type  of  boundary  conditions  will  be  classified  (References  51,  52  and 
53). 

a)  Dirichlet,  in  which  the  value  of  the  function  is  specified;  u=a. 

b)  Neumann,  in  which  the  normal  gradient  of  the  function  is  specified. 

U  =  b 

y 

c)  Mixed  (Robbins),  which  is  a  combination  of  the  above  two  types. 

U  +  bUy  =  c 

The  behavior  of  the  solution  depends  upon  the  type  of  boundary  con¬ 
ditions  applied. 

Next  we  consider  the  location  of  the  boundary  conditions,  i.e.,  either 
on  a  body  surface  boundary  or  at  a  far  field  boundary  sufficiently  removed 
from  the  body  under  investigation.  The  former  is  normally  well  defined 
geometrical ly  while  the  later  possesses  some  arbitrary  features  which  must 
be  defined  by  the  computational  aerodynamicist .  The  body  and  far  field 
boundary  conditions  will  now  be  discussed  separately. 

1 .  SURFACE  BOUNDARY  CONDITIONS 

On  the  surface  the  velocity  and  temperature  must  be  prescribed..  For 
a  viscous  fluid,  a  no  slip  condition  is  appropriate  with  no  flow  through 
the  surface.  (Although  small  amounts  of  bleed  or  suction  can  be  considered 
readi ly ) . 
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V(  77=0  =  0 

or  u,v,w(t7  =  0)  =  0 

Normally,  on  the  surface  either  a  prescribed  wall  temperature  or 
prescribed  temperature  gradient  is  used. 

T(t7=0)  =  Tw 

or  (77  =  0  =  0  Adiabatic 

<777 

More  complex  relationships  are  possible  through  consideration  of  the 
heat  transfer  process. 

Frequently  the  wall  temperature  will  not  be  known  so  that  an  estimate 
must  be  made.  Consider  the  heat  energy  balance  within  a  small  surface 
el ement . 


^  cond 

The  net  heat  into  the  element  is  obtained  through  consideration  of 
convection,  conduction  and  radiation.  This  net  heat  input  will  increase 
the  thermal  energy  of  the  element.  (Note,  under  some  ci rcumstances  , 
additional  forms  of  heat  energy  must  also  be  considered,  i.e.,  ablation, 
evaporation,  sublimation,  change  of  phase,  etc.) 
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Heat  Balance 

(^.conv- cond -*irad)dA  =  CmT  dm 

°r  •  •  -44 

h(Taw-Tw)-  Km  (Tw-Ti)-<<x  (Tw  -Tr  )  =  /?mAbCmTw 

b 

All  forms  of  heat  exchange  conceptually  can  be  grouped  into  the 
following  form: 

h  (Ta-tw)=/JmbCmTw 
or 

tw  _  _h _  _  | 

Ta'Tw  ’  '°mbCm  ~  T 

where  t  is  the  time  constant  for  the  element  to  attain  equilibrium 
temperature . 

For  convection  dominated  problems 

r . b  m  and  Ta=Taw 
P  vcp 

For  a  thin  skin  steel  model  in  a  supersonic  tunnel,  r  is  about  a 
minute.  Therefore,  adiabatic  wall  temperature  will  be  attained  in  a 
continuous  flow  tunnel.  However,  for  an  impulse  tunnel  with  running 
times  in  the  millisecond  range,  room  temperature  is  the  appropriate  value 
for  the  wal 1 . 

For  flight  application  above  M=3  the  radiation  energy  exchange  becomes 
important  and  a  "radiation  equilibrium  temperature"  is  attained. 

q  :(i  ,os  T  —  0 

tconv  trad 

4 

h  (Taw-Tw)  =  «r  Tw 

or  Tw=  — - 

!  +  €<rTw^/h 

For  this  case  Tw  is  less  than  Taw. 
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Estimates,  similar  to  these,  are  generally  acceptable  for  the  deter¬ 
mination  of  the  wall  temperature  boundary  conditions  to  be  used  in  solving 
the  Navier-Stokes  equations.  Fortunately,  the  pressure  coefficient,,  skin 
friction  and  heat  transfer  coefficients  are  not  extremely  sensitive  to  the 
value  of  the  wall  temperature.  Listed  below  are  the  dimensionless  re¬ 
lationships  for  a  laminar  boundary  layer  under  zero  pressure  gradient  for 
different  wall  temperature  (Reference  54). 


Tw/Taw 

2Re 

2St/V2Re 

tr 

0 

1.0 

.46960 

.46960 

2.591 

0.8 

ll 

II 

2.073 

0.5 

II 

ll 

1  .555 

0.4 

II 

ll 

1  .036 

0.2 

II 

ll 

0.518 

0 

ll 

II 

0 

One  additional  point  must  be  addressed  concerning  surface  boundary 
conditions.  Although  no  wall  boundary  conditions  on  either  p  or  p  are 
required  for  solving  the  differential  equations,  values  for  both  are 
needed  on  the  surface  for  the  numerical  central  difference  scheme. 

To  accomplish  this,  one  of  the  governing  equations  evaluated  at  the 
wall  may  be  utilized.  The  normal  momentum  equation  is  selected  for  this 
purpose.  The  resulting  condition  is  called  a  "compatibility  relationship" 
and  should  not  be  called  a  boundary  condition  (although  frequently  mis¬ 
labeled  in  the  literature). 

a 

p(V,  +  yVx+VVy)-y,0+<Py-rx)>.0.0 
or  Py  =  Tx 
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However,  since  F  »  t  and  |-  >>  a  simple  compatibility  condition 

<3  y  o  x 

results,  accuarate  to  order  Re-^ ,  for  the  determination  of  wall  pressure. 


Density  is  then  obtained  from  the  equation  of  state. 

Pw 

'w-  R  Tw 

This  completes  the  description  of  surface  boundary  conditions. 


2.  FAR  FIELD  BOUNDARY  CONDITIONS 

As  stated  previously,  the  specification  of  the  far  field  boundary 
conditions  depends  upon  whether  the  flow  is  subsonic  or  supersonic.  This 
section  is  intended  to  provide  some  insight  into  the  description  of  the 
outer  boundary  conditions  for  the  different  flow  classes  (Reference  55). 


To  establish  this  insight,  let  us  first  consider  the  boundary  condi¬ 
tions  for  the  simple  wave  equation. 


du  .  du 

- - rC  — 

dt  dx 


=0 


The  general  solution  of  this  equation  is  u  =  u(x-ct) 

In  the  wave  diagram  for  this  flow,  signals  are  propagated  at  speed  c 

dx 

and  therefore  boundary  condition  information  is  also  propagated  at  -rr  =  c. 


Ut  +cUx  =0 
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It  is  clear  that  in  this  case,  only  one  boundary  condition  on  x  can 
be  applied  and  must  be  applied  at  the  upstream  end  to  properly  pose  the 
problem,  i.e.,  u(o,  t) ;  Initial  conditions  are  also  required;  u  (x,o). 

It  is  also  clear  that  if  c  <  0  the  waves  propagate  in  the  opposite 
direction. 


Uf-|C|Ux=0 


0  x  L 

Again  in  this  case,  only  one  boundary  condition  in  x  is  required  but 
must  be  applied  at  the  downstream  end,  i.e.,  u  (L,  t). 

Next  consider  a  two  equation  system. 

Ut+cUx=0 

VaV° 

For  this  system,  it  is  obvious  that  the  boundary  condition  on  u  be 
placed  at  x=0,  while  the  boundary  condition  on  v  should  be  placed  at 
x  =  L.  A  set  of  initial  conditions  must  also  be  provided  for  both  u  and 
v.  In  this  form  u  and  v  are  called  the  characteristi c  variables  and 
propagated  at  speed  c  and  -a  respecti vely . 

If  we  could  find  the  characteristic  variables  for  the  fluid  dynamic 
equations,  this  would  provide  necessary  information  to  help  describe  the 
boundary  conditions. 
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To  accomplish  this,  let's  consider  the  one-dimensional  Euler  equations. 
For  most  problems  the  viscous  effects  are  minor  near  the  outer  boundary 
and  therefore,  the  inviscid  equations  are  appropriate. 


p 

pu 

pu 

+ 

/jU2+p 

pe 

t 

/>uH 

2  2 

where  H=e+P//P=^rj  +^- 


These  equations  can  be  rearranged  into  the  following  form: 

Dtt  +  u  =0 

Dt 

Du  +  a^7r  =0 
Dt 

_D  (tt-R)  =  0 
Dt 

where  7r=J_  In  p 
7 

R=  In  p 
a2=yp 

P 

The  equations  are  exact  but  non-linear.  To  obtain  the  characteri stic 
variables  we  must  assume  small  perturbations  and  linearize  the  set.  Again, 
near  the  far  field  boundary  this  assumption  does  not  lead  to  a  serious 
restri ction . 

Mow  1 et 

<^=  u  +77-  a  | 

r=  U  ~  7T  Q  | 
s  =  tt-R  (entropy) 
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Therefore,  the  inverse  relationship  are  the  following: 


R-- s  +7 r 

q+r 

U  s  1 - 

2 

In  the  new  variables,  the  three  linearized  equations  become 

+  ( u)“"  q |)  rx=  0 

st+(u,)v° 

The  characteristic  variables  are  then  q,  r  and  s  which  propagate 
at  speeds  +  a-|  ,  u^  -  a-j  and  u-j  respectively . 

It  is  also  apparent  that  q  and  s  will  always  have  positive  propagation 
speed,  (u.j  >  0  conventionally  defined  in  freestream  direction)  and  therefore 
must  be  prescribed  at  the  upstream  boundary. 

The  variable  r  possesses  either  positive  or  negative  propagation 
speed  depending  upon  whether  the  flow  is  supersonic  or  subsonic,  respectively 

Therefore,  r  is  prescribed  at  the  upstream  boundary  if  -  >  1  and  at 

ct 

the  downstream  boundary  if  —  <  1. 

d 

These  three  boundary  conditions  are  appropriate  for  one-dimensional 
unsteady  flow  and  no  additional  ones  are  needed  or  permitted.  The  numerical 
system  generally  needs  additional  information  on  the  boundaries  due  to  the 
manner  in  which  the  finite  difference  operators  are  constructed.  In  this 
case  "compatibility  conditions"  are  used  to  resolve  this  predicament  and 
should  not  be  confused  with  valid  boundary  conditions.  The  "compatibility 
conditions"  add  no  new  information  but  merely  redescribe  the  available 
information,  e.g.,  insuring  that  the  governing  equations  are  satisfied  at 
the  boundary.  For  the  one-dimensional  Euler  equations,  we  w-ill  merely 
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reuse  the  characteristic  equations  at  the  boundary  to  resolve  the  undefined 
variables.  Below  is  a  summary  of  the  appropriate  boundary  conditions  for 
either  supersonic  or  subsonic  flow. 

3.  SUMMARY  ONE-DIMENSIONAL  FLOW 
Far  Field  Boundary  Conditions 
Subsonic 


*2  L 

v*. 

ht+<u+al,jJ2=0 

[rj+fu-olrJ^O 

r=r2 

s=s, 

fst+(u)s  ]  =0 

1 

L  T  ^2 

Supersonic 

x  =0 

x0  =  L 

[it+(u+o,<ix]2=0 

r=ri 

[rt  +  (u-(i)r,]2»0 

S  =  S| 

[st+(u)s  ]  =0 

L  *  ^2 

The  following  figures  show  the  implementation  of  these  boundary 
conditions  for  the  Navier-Stokes  equations  for  MOT  =  0.5  uniform  flow 
state . 

Initially  a  square  wave  is  inserted  into  the  middle  of  the  flow  field 
and  Figures  8  and  9  show  the  propagation  of  the  character!' stic  variables 
(q,  r,s)  at  different  time  levels.  These  are  the  correct  boundary  con¬ 
ditions  in  this  case  and  no  reflections  at  the  boundary  occur.  In  Figure  10 
a  case  is  displayed  for  the  inappropriate  boundary  conditions  (but  commonly 
employed)  of  p  ,  U  ,  P  given  upstream  and  p  Uv ,  P  downstream.  Note  in 

00  OO  CO  X  A  X 
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Wave  Diagram  for  Initial  Condition  Disturbances 
(Correct  Boundary  Conditions) 
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this  case  that  waves  continue  to  reflect  from  the  boundaries  resulting  in 
1 arge  errors . 

4.  BRANCH  CUT  BOUNDARY  CONDITIONS 

Frequently  within  the  computational  domain  a  branch  cut  is  inserted 
which  joins  the  inner  and  outer  boundaries  dividing  the  flow  field  into 
separate  sections.  For  example,  for  symmetrical  configurations  a  branch 
cut  is  located  on  the  plane  of  symmetry  which  permits  one  to  solve  only 
one-half  of  the  problem,  thereby  saving  one-half  of  the  computer  resources. 
Along  these  branch  cuts  boundary  conditions  must  be  applied.  Two  types 
occur,  i.e.,  symmetric  and  periodic. 

5.  SYMMETRY  PLANE 

For  this  situation,  no  flow  is  permitted  through  the  plane  and  all 
gradients  normal  to  the  plane  must  vanish 

v  =  0  and  dLT  -q 

dn 


6.  PERIODIC  CONDITIONS 

For  the  case  of  an  arbitrary  artificial  boundary  located  in  the  field, 
for  example,  encountered  in  a  cascade  of  turbine  blades,  all  properties 
must  be  continued  across  the  cut.  Due  to  the  periodic  nature  of  this 
situation  the  following  boundary  condition  is  appropriate 


Where  stations  1  and  N  in  the  transformed  plane  are  geometrically 
coincident  stations  in  the  physical  plane  representing  the  multivalued 
features  of  the  conf i guration . 

7.  CLASSIFICATION  OF  PARTIAL  DIFFERENTIAL  EQUATIONS 

Non-linear  partial  differential  equations  can  be  classified  according 
to  the  type  of  subsidiary  condition  that  must  be  imposed  to  give  a  well -posed 
problem.  Consider  a  non-linear,  second  order  quasi -linear  partial  differential 
equation. 
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2  2  2 

A  ^4*  +  2B  j-j’  +  C-^-  =  <i>(ii  ,ux,uy,x,y) 
2  dxdy  •>  2 


dx 


dy 


where 


A,B,C  =  functions  of  x,  y  only 


Also  assumed  valid  in  the  x,y  plane  are  the  total  derivative  definitions. 

dux5^-  (ux)dx+r-(ux)dy 
dx  dy 

dijy=^(uy)dx  +  “(uy)dy 

This  constitutes  three  equations  for  the  determination  of  u  ,  u  and  u 

^  xx  xy  yy 


Determinant5  D= 


A  2B  C 

dx  dy  o 

0  dx  dy 


Two  families  of  characteristics  (real  or  complex  conjugates)  curves  exist 
on  which  0=0.  This  relation  is  known  as  the  equation  of  characteristics . 


Ady2-  2Bdxdy+Cdx2=0 


or 


dy 


dx 


B± 


2 

B  -AC 


Three  classes  of  equations  have  been  identified  dependent  upon  the  con¬ 
dition  of  the  radical  (Reference  5 

1 .  Elliptic  32  -  AC  <  0 

The  characteri sti cs  are  complex  conjugates 

2 .  Parabol ic  B2  -  AC  =  0 

Only  one  real  family  of  characteri sties  exist. 
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3.  Hyperbolic  a  -  AC  >  0 
Characteristic  are  real 

Examples  of  the  three  classes  are  as  follows: 

1.  Laplace  Equation  for  two-dimensional  flow. 

/<£  dZ$ 

—  + — -  =0 

ax2  ay2 


A=  I 


B  =  0  B  -AC  =  -|<  0  Elliptical 
dy 

c=l  i=±i 


2.  Heat  conduction  in  thermodynamics 

dt  2 

dx 


A=0 


B=0  B  -AC  =  0  Parabolic 

C’“  |2.0 

3.  Vibrating  string  in  mechanics. 

2  2 
<3_y  _  2  £  y 

2  °  2  "  U 

at  ax^ 


As  I 

2  2 

B=0  B  -AC  =  a  >0  Hyperbolic 

C=-a2,±x  =  +g 
dt 
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The  description  of  the  outer  boundary  conditions  for  flow  fields  is  quite 
different  depending  upon  the  classification  of  the  type  of  partial 
differential  equations. 
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SECTION  XV 

GRID  GENERATION  PROCEDURE 

An  automated  procedure  is  needed  to  generate  a  body-oriented 
coordinate  system  for  arbitrary  geometries.  Three  types  shall  be  dis¬ 
cussed,  i.e.,  algebraic,  elliptic  and  hyperbolic;  which  are  typical  of 
the  types  presently  in  use.  This  is  an  area,  however,  where  improvements 
are  required  in  order  to  upgrade  the  overall  efficiency  of  computational 
aerodynamics . 

1 .  ALGEBRAIC  METHOD 

The  homotropy  method  of  R.  Smith  (Reference  55)  is  an  example  of  an 
algebraic  method.  In  this  approach,  geometric-constructs  are  used  to 
define  a  grid  without  resort  to  differential  equations.  First,  define 
the  body  contour  with  grid  points  located  at  constant  increments  of  arc 
length  (Figure  11).  Label  this  curve  n=0.  Next,  construct  an  outer 
boundary  with  the  same  number  of  grid  points  as  on  the  body  and  also  at 
constant  increments  of  arc  length.  Label  this  curve  n=N.  These  two 
curves  need  not  be  of  equal  arc  length  but  they  must  possess  the  same 
number  of  grid  points.  Now  connect  the  first  point  of  the  n=0  curve 
with  a  straight  line  to  the  first  point  of  the  n=N  curve.  Next, 
sequentially  connect  all  the  points  between  the  two  curves  with  straight 
lines.  Label  these  straight  lines  5  =  0,  1,  2  ...M  in  numerical  sequence. 
Now  divide  all  g  lines  into  N  steps  of  similar  proportion.  Any  proportion 
is  feasible  although  a  systematic  regular  variation  produces  the  best 
results.  (This  maintains  better  behavior  of  the  higher  derivatives  of 
the  metrics).  Now  connect  these  points  to  form  the  family  of  n  =  constant 
curves.  From  this  constructed  network  a  one-to-one  relationship  of 
x(£,n)  and  y(c,n)  can  be  determined  for  which  the  metrics  can  be  computed 
numerically.  An  example  of  this  grid  is  shown  in  Figure  12. 
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2.  ELLIPTIC  METHOD 

The  pioneering  method  of  J.  Thompson  (Reference  57)  in  the  use  of  a 
general  transformation  procedure  was  responsible  for  early  successes  in 
the  field.  This  procedure  (Reference  46)  also  involves  establishing  an 
inner  body  contour  (n=N).  On  these  contours  a  grouping  of  (x,y)  points 
is  selected  to  define  a  closed  contour;  C-j  (x,y)  on  the  inner  and 
C^(x,y)  on  the  outer  contour  (Figure  13) 

A  simple  mapping  is  desired  to  determine  the  interior  points  with 
constraints  that  a  minimum  or  maximum  not  occur  in  the  interior  (to 
maintain  single-valuedness)  and  also  that  coordinate  lines  of  the  same 
family  not  intersect.  The  elliptic  Laplace  equation  contains  these 
desired  features. 

_2  _ 

V  77=03'17  +  T) 

'  'xx  '  yy 

This  equation  with  boundary  conditions  on  and  completely 
define  the  problem. 

A  second  family  is  also  generated  by  the  Laplacian  with  periodic 
boundary  conditions  at  an  arbitrary  branch  cut. 

V2£= 0 

To  solve  these  equations  we  resort  to  the  inverse  transformation 
and  exchange  dependent  and  independent  variables.  It  will  be  seen  that 
numerical  solution  in  the  transformed  plane  is  more  convenient. 

The  chain  rule  states  that 

±s,  i_.  J__ 

dx  x  <3£  ^x  dr, 

—  —  +T>  — 

ay  y  d£  1y  dr, 
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where 


£  .j“'v  £  =-J  *x 

%x  yv  v 

77  ^  V  77  *  J  X ^ 

7x  7£  'y  £ 


J  =  x. y  -  x_y- 

£,77  ^  £ 


Applying  the  inverse  transformation  to  the  Laplacian  equation  pro¬ 
duces  the  following: 

2  2 

J  V  £  =  (y  y  *-y  y  +  x  x  .-x  x  ) 

V  7?£  ^7777  77  77£  £  7777 


■aj  *J^ 


and 


,2„2 


J_V>=(-yi7y^  +  yjy^-x,Xj£  +  x?xfi))+/3j‘lJf-yj'lJ 


a  =  xf;+yf; 
77  '77 


where 


/?=  xA  x  +y„  y 
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Combining  these  relationships  produces  the  transformed  Laplacian  for  x 
and  y. 

-J2(x^V2C+x7?V2'7)=0=ax^-2/3  x^+yx^ 

-  J  2(y$V2£  +  y^V  %  )=0=ay^-  Z&y^+yy^ 

These  two  sets  of  coupled  elliptic  partial  differential  equations  require 
numerical  solution.  Two  methods  have  been  used.,  i.e.  SOR  (successive-over 
relation)  and  ADI  (alternating  direction  implicit  method). 

3.  SOR  SOLVER 

The  method  of  successive-over-relaxation  (SOR)  was  developed  by  Young 
(1954).  It  is  an  iteration  method  to  relax  the  equation  from  some  initial 
guess  by  driving  the  error  residuals  to  zero  at  each  point.  The  term 
"over-relaxing"  implies  applying  a  larger  correction  than  the  standard 
relaxation  calculation  produces  in  order  to  accelerate  convergence. 

To  demonstrate  the  procedure  consider  the  original  Laplace  equation. 

V2^  0 

In  finite  differences  it  becomes 

A+i,r2A,j+£i-i,j  ,  ci,i+r2A,j+A,i-i 

Ax2  Ax2 

If  Ax  =  Ay  for  simplicity,  then  at  any  iteration  cycle  the  residual 
becomes 
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These  residuals  are  calculated  at  each  point  in  the  field.  Corrections 
are  then  applied  in  a  systematic  fashion  for  the  next  iteration.  Various 
methods  have  been  developed  to  improve  convergence.  SOR  uses  the  following: 


i.l  4 


where  n  =  iteration  level 

l<  u»  <  2=  relaxation  factor 


The  method  is  continued  until  the  residuals  are  driven  to  sufficiently 
small  valves. 

4.  ADI  SOLVER 

The  alternating  direction  implicit  method  (ADI)  was  introduced  by 
Peaceman  and  Rachford  (1955)  (Reference  31).  The  method  splits  the  equations 
into  two  one-dimensional  parts  and  uses  the  efficient  tridiagonal  solver 
in  .each  direction  al ternatively  (References  58  and  59).  First,  sweep  in 
one  direction  while  holding  the  derivatives  constant  in  the  other  and 
then  reverse  the  procedure  to  complete  the  cycle. 


Step  1 


C/Z-2in+'/ZHn+!/2-A^i^ 

i  +  l,J  i,J  i-l,J  ay2 


*J 


i-l,j 


Step 


2  C-<l+Cr-A>2(ii,n+l'2 

U  +  l  M  l»J~l 


Again,  the  alternative  sweeps  are  continued  until  convergence  is  achieved 


5.  HYPERBOLIC  METHOD 

In  external  aerodynamics,  the  location  of  the  outer  boundary  need 
not  be  specified;  it  only  need  be  far  removed  from  the  inner  boundary. 
Hyperbolic  methods  (Reference  5)  may  be  used  in  this  case  as  developed 
by  J.  Steger  (1980).  We  seek  a  grid  composed  of  constant  ?  and  n  lines 
given  initial  data  along  n  =  0  on  the  body  contour.  A  set  of  partial 
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differential  equations  are  sought  to  generate  a  smoothly  varying  mesh 
such  that  grid  lines  of  the  same  family  do  not  intersect  or  coalese, 

These  equations  may  be  obtained  from  two  conditions,  i.e.  an  orthogonal i ty 
condition  and  a  geometric  constraint.. 


WW° 


'J<v> 


(Cauchy-Riemann) 
Orthogonal ity  Condition 

=  Area  Constraint 


Cramer's  Rule  may  be  used  to  solve  for  x  and  y 
J  n  n 


-Jy, 


J  ^ 

J  x. 


V  V  2 


xl+ye 


■ye 


A 

x|  +  y 


2 

t 


These  equations  are  hyperbolic  and  can  be  marched  in  the  n  direction. 


The  advantages  of  the  method  is  that  it  is  fast,  orthogonal,  auto¬ 
mated  and  clustering  can  be  controlled  by  varying  J.  It  can  only  be  used, 
however,  when  the  outer  boundary  need  not  be  specified  (Figure  14). 
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Figure  13.  Sketch  of  Physical  and  Computational  Plane 
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SECTION  XV 

FLUID  DYNAMIC  STABILITY 

In  a  previous  section  we  analyzed  the  stability  of  the  numerical 
method.  Criteria  were  established  to  insure  that  a  stable  algorithm  was 
utilized  in  solving  the  fluid  dynamic  equations.  The  physical  flows 
however  can  also  exhibit  an  instability  due  to  natural  causes.  It  is  the 
purpose  of  this  section  to  identify  the  situation  under  which  real 
instabilities  can  exist  in  order  to  help  discern  the  difference  from  an 
unphysi cal  numerical  .instability. 

To  demonstrate  the  procedure  a  common  fluid  dynamic  instability  will 
be  investigated  entitled  the  "Rayleigh  Instability"  (References  26  and  60). 

Examine  the  incompressible,  inviscid  (Euler)  equations. 

ux+vy=0 

Ut+iiUX+VUy-PX^ 

vt.+u \+vy-py/p 

These  last  two  equations  may  be  combined  to  eliminate  pressure  by 
introducing  vorticity. 

NOTE:  to  f  0 

Rotational  flow  is  considered  here. 

VxV  =  (V  -U  )k 

—  —  —  x  y  - 

1 .  PARALLEL  FLOW 

Assume  that  the  flow  can  be  represented  as  a  disturbance  from  a 
steady-state  parallel  shear  flow  as  follows: 

U=U(y)+U'(x,y,t) 

V=0+V'(x,ytt) 


D<£  =0 

"d? 

where 
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Hence  .  .  _  _ 

Vx-Uy-(Vx-Uy)-Uyaw'-Uy 


The  governing  equations  become 


V-V=U  x+Vy=0 

7^  =  0,' +  Uo/ -U  V'+(H.0.T.)  =  0 
Dt  t  x  yy 

The  boundary  conditions  are  that  and  V"  vanish  at  ±  ». 

The  governing  equations  are  linear  and  possess  the  following  solution 

U/=U(y)eia(x'C,) 

V'^(y)eia(x-C,) 

where 

A 

U,<£  and  c  are  complex 
a  =  real  (wave  number) 

Therefore 

iaU  +  <£y=0 

<u'-(ia*-Uy)eia(x‘c,) 

(U-C)(iaUy+a2*)  +  Uyy«#.=0 

Eliminating  U  produces  the  Rayleigh  Equation. 

fyy-‘a2+fe>f-0 

with  boundary  condition 


<£(±CO)  =  0 
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Given  U  and  a  these  boundary  conditions  can  only  be  achieved  for  specific 
values  of  C  =  Cr  +  id.  This  eigenvalue  equation  was  first  studied  by 
Rayleigh  in  1880.  He  derived  a  simple  criterion  for  the  condition  under- 
which  the  flow  was  unstable,  i.e.  d  >  0 

2.  CONDITION  FOR  INSTABILITY 

Multiply  the  Rayleigh,  equation  by  the  conjugate  of  <f>  (denoted  by 
(j>*)  and  integrate  over  the  entire  domain. 

_ £  [$"<£* ~  a 2  4> 4>*~ </><£*]  dy*0 


Manipulating  the  terms  to  obtain  the  real  and  imaginary  parts  produce  the 
following. 


^/ [(<£'<£*) 


UlO-cW, 

(U-C)  (U-C*)-* 


Since  the  first  term  vanishes  due  to  boundary  conditions 

=0 

and  the  second  and  third  terms  are  real,  only  the  last  term  contains  any 
imaginary  part. 


Hence 


i  C 


cp  0 "<p  <£* 

]  f  —  2  2 

-co  (U-Crr+Cf 


dy  =  Q 


This  relationship  can  be  satisfied  in  two  ways,  i.e.  either  d  =  0  or 
U"  =  0  somewhere.  The  later  is  the  essential  condition  for  an  instability. 
This  implies  that  the  velocity  profile  exhibits  an  inflection  point.  This 
deduction  is  entitled  Rayleigh's  Second  Theorem. 
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By  inserting  a  value  for  U  =  U(y)  into  the  Rayleigh  equation,  eigen¬ 
values  for  C  =  C(«)  can  be  determined.  This  was  accomplished  by  Verma, 
Hankey  and  Scherr  (Reference  61)  for  the  series  of  separated  boundary 
layer  profiles  obtained  from  the  Lower  Branch  solution  of  the  Falkner-Skan 
equations.  A  plot  of  these  results  indicate  that  all  the  velocity  profiles 
with  inflection  points  are  unstable  (as  indicated  by  Rayleigh's  second 
theorem),  however,  only  for  a  small  range  of  frequency  «<s  =  .  In 

addition,  large  values  of  Cl  are  evident  indicating  a  severe  instability 
will  occur. 
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SECTION  XVII 
TURBULENCE  MODELS 


Although  turbulence  is  a  self-excited  oscillation  and  predictable 
using  the  Navier-Stokes  equations,  the  scale  of  the  smallest  eddy  would 
require  extremely  small  grid  sizes,  thereby  rendering  the  computation 
impractical.  As  a  consequence,  turbulence  is  treated  as  a  fluid  property 
and  empirically  added  to  the  equations.  In  this  section  we  will  develop 
this  concept. 


Beginning  with  the  two-dimensional,  unsteady,  incompressible  Navier- 
Stokes  equations  we  shall  derive  the  Reynol ds-averaged  equations  (Reference 
26): 


WFy=° 


0 

u 

V 

u= 

U 

;E= 

2 

U  —  O' 

\\/p 

i  F« 

UV-T//J 

V 

UV-T/yO 

v2"°22//j 

These  variables  are  considered  to  be  fluctuating  in  time  about  a 
well  defined  mean. 


where 


u=u  (x,y)+u' (x^,!) 

v=v(x,y)  +  v'(x,y,t) 
p=p(x,y)  +  p/(x,y,t) 


t+P 


udt 


and 


♦+P 

f‘  Vdt=o 
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where  the  period  over  which  the  average  is  accomplished  is  large  compared 
to  the  period  for  the  lowest  frequency  component  of  concern. 

p«  f  min  (~  10  Hz  for  example) 

Integrating  and  determining  the  mean  produces  the  following: 

I  ,t+P 

p,  aJ*  +  Ex+Fy)dt=0 

The  linear  terms  (u,  v,  p,  t...,  a^.),  merely  produce  the  mean  values. 

The  non-linear  terms  produce  additional  terms. 


I  rt+P i,2 I  ft+P-2.  /.  >2,  ..  -2  ,2 

p/^  u  dt  =  p /f  (u  +2uu+u  )dt=u  +u 

similarly 

I  /-t+p  2  -2  ~~/2 

=  /  v  dt=  v  +  v 

P  t 

I  t+P  —  — 

p /  u  vdt  =  u  v  +  (/  v/ 


Therefore 


0 

U 

V 

u« 

IT 

;  E= 

-2  ^2  — — 

U  +  U  —  cr  II//J 

i  F  = 

U~  +  u/v/-  a/p 

V 

u  v  +  u/  v'-f/p 

V2  +  V  Z~cr22/p 

The  additional  terms  all  appear  next  to  one  of  the  stress  terms.  These 
new  terms  are  entitled  "apparent  stresses"  or  Reynolds  stresses.  It  is 
therefore  convenient  to  redefine  these  stresses  for  turbulent  flow  to 
make  the  equations  identical  in  form  to  the  laminar  equations. 


Xrb=_S+W'-+2^_x'',U 

rlo  =^(uy+vx)— /ju'  v7 
‘^turb  _ 

r22turb  -P+XVT  +  2^y-/’v/‘ 
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By  analogy  with  the  description  of  viscous  stresses  created  by  the 
molecular  viscosity,  we  can  define  a  turbulent  (or  eddy)  viscosity  as 
follows : 

~P  |^(ITy  +7x) 

-P  u/2s«n(-2/3V-v+2ux) 

-/»v  ,2s«22(-2/3V'7+2  vy) 


It  can  be  shown  that  e  is  always  positive  defined  in  this  manner  to 
prevent  violation  of  the  second  law  of  thermodynamics.  Since  insufficient 
empirical  information  is  available  to  evaluate  these  different  eddy 
viscosity  coefficients  they  are  equated  to  each  other.  The  largest  term 
of  engineering  importance  is  U  and  requires  the  greatest  attention. 
The  remaining  terms  generally  contribute  little  and  need  not  be  evaluated 
accurately. 

Therefore  *12*  *11*  *22 

The  magnitude  of  the  apparent  stresses  overwhelms  the  molecular 
viscosity  terms  (except  on  the  surface)  since 

100  to  1000  for  I06<  Re  <  I08 

Therefore,  the  Reynolds-averaged  equations  are  obtained  simply  by 
replacing  y  by  e  in  all  stress  terms. 


For  compressible  flow  we  include  the  density  fluctuations  and  the 
energy  equation,  i.e.,  p=p+p' 

However,  it  can  be  shown  that  the  p "  variation  produces  little  change 
in  the  equations  up  to  M=5.  This  is  due  in  part  to  the  fact  that  the  U" 
disturbance  is  still  subsonic  up  to  M  ~  5.  In  the  energy  equation  the 
thermal  conductivity  is  replaced  by  the  turbulent  conductivity. 


K  =  Kf 
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which  is  evaluated  by  means  of  the  eddy  viscosity, 
Prandtl  number. 


lam 


A-Cp 

K 


0.72 


£• ,  and  turbulent 


_  €  Cp 

Pr+  u=l r*  3  1.0 
turb  Kt 


Therefore,  the  compressible,  turbulent  Reynol ds-averaged  equations  are 
identical  to  the  compressible  laminar  Navier-Stokes  equations  with  new 
values  for  the  transport  properties. 


1.  EVALUATION  OF  EDDY  VISCOSITY 

Thus,  there  is  no  difference  between  calculating  laminar  or  turbulent 
flow  except  for  the  calculation  of  e.  In  this  section,  we  will  examine 
simple  turbulence  models  (Reference  62).  Although  complex  relationships 
involving  a  system  of  partial  differential  equations  with  27  unknown 
coefficient  have  been  derived  to  evaluate  the  Reynol ds-stresses ,  they  have 
not  produced  as  originally  advertised.  The  present  method  in  vogue  today 
is  to  use  simple  algebraic  models  and  adjust  the  constants  for  the  special 
cases  under  consideration.  Hopefully,  a  pattern  will  evolve  as  more 
experimental  data  and  numerical  comparisons  are  produced. 


2.  BOUSSINESQ  MODEL 


By  forming  a  dimensionless  turbulent,  Reynold's  number  for  turbulent 
boundary  layers  a  correlation  was  determined. 


Pu  S 

Returb=  -T- 


=  60 


3.  CLAUSER  LAW  OF  THE  WAKE 

It  was  later  found  that  using  the  incompressible  displacement  thick- 

k 

ness,  6.,  for  the  length  scale,  the  outer  80S  wake-like  region  of  com¬ 
pressible  turbulent  boundary  layers  could  be  correlated. 
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€outer=  ^2^^® 


where  Kg^SO)-1  =.0168 


4.  VON-KARMAN  LAW  OF  THE  WALL 

In  the  inner  20%  e  diminishes  from  the  above  value  due  to  the 
presence  of  the  wal  1  . 


Experiments  indicated  a  logarithmic  velocity  variation  near  the 
wall  from  which  the  eddy  viscosity  could  be  deduced. 


«a/»(K.y) 


du 

dy 


where  K 


0.4 


The  shear  stress  is  nearly  constant  in  the  vicinity  of  the  wall  for 
zero  longitudinal  pressure  gradient  (flat  plate). 


dr 

dy 


dp 

dx 


Therefore  r  =  r  =  « 

w 


du 

dy 


»7,(v$!)2-u*z 


Integrating  this  expression  reaffirms  the  logarithmic  velocity 
profile  variation. 


JJ  1 
U*'K, 


Experiments  show  that  C=5  for  smooth  flat  plates. 
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5.  VAN  DRIEST  DAMPING  FACTOR 

Refined  measurements  near  the  wall  identified  the  existence  of  a 
laminar  sublayer.  Van  Driest  developed  an  expression  to  join  the  laminar 
region  to  the  law  of  the  wall  region. 


The  following  damping  factor,  D,  was  inserted  into  the  length  scale. 

-yU* 

_  .  i/A 
D= 1  — e 

where  A=26 

The  inner  value  for  the  eddy  viscosity  becomes  the  following: 


inner 


=  />(K|yD)‘ 


du 

dy 


6.  CEBECI-SMITH  MODEL 

A  combination  of  the  above  relations  is  called  the  Cebesci -Smi th 
Model.  A  two-layer  turbulence  model  is  used  for  the  inner  and  outer 
eddy  viscosity.  Both  are  programmed  and  computed  separately,  however, 
the  lesser  value  of  the  two  is  used  in  each  region. 

7.  BALDWIN-LOMAX  MODEL 

Since  the  two-layer  model  is  somewhat  awkward  in  joining  two  separate 
functions,  a  unified  method  is  preferred.  In  addition,  relating  the 
turbulence  to  vorticity  is  regarded  as  fundamental  by  some  investigators. 
The  Bal dwi n-Lomax  model  accomplishes  these  features. 

2 

€  =  /0(DI)  I  oil 

where  1=  .088  tan  h  (4U-) 

.08  s 


8.  FAR  WAKE  MODEL 

For  far  wakes  the  form  of  the  law  of  the  wake  is  used,  however,  the 
coefficient  increases  by  a  factor  of  4. 
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<wakes-064^UeSi 

Empirical  correlations  are.  used  to  join  s  in  the  various  regions 
of  the  flowfield. 
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SECTION  XVIII 
SUMMARY 

Included  herein  is  an  introduction  to  Computational  Aerodynamics 
in  which  the  major  topics  are  addressed.  The  purpose  of  this  report 
is  to  provide  a  foundation  for  those  just  entering  the  field.  It  is 
intended  that  additional  sections  be  added  in  the  future  as  further 
developments  evolve  in  the  subject  area. 
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